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ABSTRACT 

• 

- _ -f  Large  sets  of  coupled,  nonlinear  equations  arise  in  a  number  of  disci¬ 

plines  in  connection  with  computer  based  models  of  physical,  social  and 
economic  processes.  Solutions  for  such  large  systems  of  equations  must  be 
effected  by  means  of  digital  computers  using  appropriately  designed  codes. 

This  paper  addresses  itself  to  the  critically  important  problem  of  how 
sensitive  the  solutions  are  to  variations  of,  or  inherent  uncertainties  in, 
the  parameters  of  the  equation  set.  We  review  here,  and  also  present  further 
developments,  of  our  statistical  method  of  sensitivity  analysis.  The  sensiti¬ 
vity  analysis  presented  here  is  nonlinear  and  thus  permits  one  to  study  the 
effects  of  large  deviations  from  the  nominal  parameter  values.  In  addition, 
since  all  parameters  are  varied  simultaneously,  one  can  explore  regions  of 
parameter  space  where  several  parameters  deviate  simultaneously  from  their 
nominal  values.  /  .  W 

0(1*  I Op* C\  »  A  ^ 

We^e\H*l^\here/ita  theory  of/oer'method  of  sensitivity  analysis, 
then  detail  the  method  of  implementation  and  finally  present  several  examples 
of  its  use  to  date.  ^ 


I .  INTRODUCTION 


Sets  of  coupled,  non-linear  equations  arise  in  a  number  of  dis¬ 
ciplines  in  connection  with  computer  based  models  of  physical,  social 
and  economic  processes.  These  sets  of  equations  may  be  differential, 
integral  or  algebraic.  They  arise  in  such  widely  different  fields  as 
reaction  kinetics,  combustion,  air  pollution,  weather  forecasting,  upper 
atmosphere  phenomena,  seismic  analysis,  operations  research,  systems 
analysis,  economics,  etc.  These  model  systems  may  contain  as  many  as 
100  equations,  a  very  large  number  of  parameters  (in  the  form  of  coupling 
coefficients  such  as  rate  coefficients,  transport  coefficients,  economic 
coefficients,  etc.)  and  a  very  large  number  of  output  variables.  Solutions 
for  such  large  systems  of  equations  must  be  effected  by  means  of  digital 
computers  using  appropriately  designed  "codes". 

As  is  well  known,  the  computer  solution  of  such  large  sets  of 
equations  can  be  quite  expensive.  Even  after  such  solutions  have  been 
achieved  one  is  still  faced  with  a  critically  important  problem:  How 
sensitive  is  the  solution  to  variations  of,  or  Inherent  uncertainties  in, 
the  parameters  of  the  equation  set?  This  problem  of  "sensitivity"  Is 
central  to  the  understanding  of  the  behavior  of  systems,  and  of  the  models 
representing  such  systems,  which  contain  a  large  set  of  coupled  equations. 
It  is  clearly  Important  to  know  how  sensitive  the  output  variables  are  to 
changes  of,  or  uncertainties  in,  the  parameters  and  which  of  the  variables 
are  sensitive  (or  not  sensitive)  to  which  of  the  parameters.  Until  this 
Information  Is  available,  any  proposed  model  must  be  suspect  as  a  valid 
representation  of  the  real  system.  Furthermore,  the  accuracy  to  which  the 
model  parameters,  i.e.,  coupling  coefficients,  need  to  be  determined 
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via  calculation  or  experiment  depends  upon  the  sensitivity  of  the  output 
variables  to  the  value  of  the  parameters.  And  finally,  any  desired 
optimization  of  various  output  variables  with  respect  to  the  coupling 
coefficients  requires  a  knowledge  of  the  sensitivity. 

Any  attempt  to  determine  the  "sensitivity"  by  solving  the  set  of 
equations  over  and  over  aga’in,  varying  one  parameter  at  a  time  over  a 
series  of  values  while  holding  all  the  other  parameters  fixed  at  some 
specific  values  becomes  prohibitive  in  time  and  expense  for  the  large 
systems  discussed  here.  This  is  readily  demonstrated  by  a  simple 
calculation.  For  a  model  system  of  many  coupled  differential  equations 
with  n  parameters  and  m  output  variables,  the  above  procedure  for  z 
different  values  for  each  of  the  parameters,  would  require  zn  Integrations 
for  each  of  the  m  variables,  i.e.,  a  total  of  m(z)n  integrations.  For  m, 
z  and  n  large,  not  only  will  the  computations  be  prohibitively  expensive, 
but  the  printouts  would  be  so  numerous  that  the  analysis  of  the  results 
themselves  will  be  a  major  problem. 

In  response  to  this  need,  we  have  developed  a  statistical  method  for 

the  sensitivity  analysis  of  large  systems  of  coupled  non-linear  equations. 

1  2 
The  theory  ,  its  application  to  several  test  kinetics  systems  and  an 

analysis  of  the  approximations^  have  already  been  presented  and  the  reader 

is  referred  to  these  publications  for  various  details  which  may  not  be 

covered  in  this  paper. 

The  purpose  of  this  article  is  to  recast  and  further  develop  our  method 

of  sensitivity  analysis  Into  a  form  which  has  a  number  of  advantages  over 

11-31 

the  previous  formulation'  .  Jhis  new  formulation  permits  us  to  discuss 
sensitivity  analysis  from  a  more  familiar  and  direct  point  of  view.  A 
further  advantage  is  that  the  relation  to  more  conventional  sensitivity 
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analysis  is  now  easily  obtained.  Finally,  this  article  is  user  oriented. 

We  present  here  all  the  steps  required  for  the  readers  to  apply  this  method 

of  sensitivity  analysis  to  their  own  systems. 

It  is  very  important  to  point  out  here  that  the  implementation  of  our 

method  of  sensitivity  analysis  is  very  simple  even  though  the  theoretical 

analysis  presented  here  appears  to  be  quite  involved.  We  will  demonstrate 

this  in  Sections  III  and  IV  where  we  discuss  the  implementation  of  the 

method  and  present  some  specific  examples. 

A  mathematical  model  of  a  physical  system  is  a  programmed  computer 

algorithm  which  returns  a  prediction  y  (y  may  be  a  vector)  for  anjr. 

physically  realizable  values  of  the  parameters  k  and  constants  c,  and 

over  any  physically  meaningful  range  of  values  of  the  independent  variables 

x.  Such  a  model  may  return  nonsense  for  certain  combinations  of  values  of 

parameters  and  independent  variables;  but  we  assume  that  the  usual  theory 

vs.  experiment  checks  have  already  been  carried  out  on  the  model,  so  that 

♦ 

gross  deficiencies  are  not  evident. 


♦Clearly,  If  these  conditions  are  not  met,  there  is  a  problem  at  a  level 
more  basic  than  that  which  calls  for  a  sensitivity  analysis.  We  might 
mention  that  our  experience  has  shown  that  computer  algorithms  frequently 
have  been  checked  out  only  for  a  limited  number  of  specific  parameters  values, 
and  not  over  the  broader  range  of  values  for  which  the  model  presumably  is 
valid.  In  order  to  apply  the  sensitivity  method  discussed  here  (or  any 
other  method  of  sensitivity  analysis),  it  is  necessary  to  "tune"  the 
algorithm  to  cover  the  broad  range  of  values  which  the  model  is  presumed 
to  represent.  In  unfavorable  circumstances,  this  may  require  more  work  than 
the  sensitivity  analysis  itself. 
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A  model  typically  is  characterized  or  controlled  by  "parameters", 
as  shown  in  figure  1.  In  the  abstract,  we  might  regard  these  parameters 
as  being  a  subset  of  the  independent  variables  except  for  a  very  important 
distinction.  The  true  independent  variabl-es  always  cover  a  range  of  values 
(possibly  infinite)  during  a  single  run  of  the  model.  The  simplest  example 
of  such  an  independent  variable  is  time  in  dynamic  processes,  but  space 
and  many  other  variables  could  be  Independent  in  specific  problems.  The 
parameters,  on  the  other  hand,  have  unique  values  during  the  course  of  a 
single  "run"  of  the  model,  although  It  may  be  necessary  to  vary  these  values 
from  one  run  to  the  next.  The  need  for  such  variation  may  arise  from  any 
of  several  possible  sources,  several  of  which  are  indicated  in  Fig.  1.  For 
example,  a  "physical"  parameter  certainly  has  an  unique  value,  but  this 
exact  value  may  not  be  known  because  of  limitations  of  Information;  only 
a  range  or  distribution  of  values  may  be  known.  As  another  example,  a  para¬ 
meter  may  be  controllable  in  a  particular  physical  circumstance  only  to 
within  some  range,  e.g.,  the  impedance  of  a  variable  element  of  an  electroni 
circuit,  or  the  mass  loading  of  a  spring  In  a  mechanical  system. 

Separate  also  in  definition  from  the  independent  variables  and  para¬ 
meter  are  the  fixed  constants  of  the  model,  which  do  not  vary  within  the 
context  of  the  class  of  problems  of  interest  to  the  model  user,  and  whose 
values  can  be  precisely  specified.  It  should  be  noted,  of  course,  that 
what  Is  a  fixed  constant  in  the  context  of  one  situation  might  be  a  para¬ 
meter  in  the  context  of  another  situation;  the  distinction  depends  on 
the  particular  case  on  hand. 

The  fact  that  the  parameters  can  take  on  a  range  of  values  suggest 
that  a  statistical  approach  to  sensitivity  analysis  is  appropriate.  Instead 


of  considering  the  effect  on  the  output  functions  of  one  at  a  time  variations 
in  each  of  the  parameters  as  in  a  "brute  force"  method,  we  will  construct 
outputs  averaged  in  one  operation  over  probability  distributions  of  all 
the  parameters.  The  distribution  of  the  parameters  can  arise  because  of 
experimental  uncertainties  or  theoretical  approximations,  because  of 
"ignorance"  of  the  value  wfthin  certain  reasonable  bounds,  or  might  represent 
upper  and  lower  limits  due  to  "stops"  on  the  physical  controls  of  the  systems 
being  modeled. 

Our  method  of  sensitivity  analysis  proceeds  by  relating  the  probability 
distribution  of  each  parameter  to  a  frequency  and  one  new  parameter  s  which, 
as  s  varies,  carries  all  the  parameters  through  their  range  of  variation. 

The  parameter  s  is  varied,  and  the  frequencies  chosen  in  such  a  way  that 
the  output  variables  at  any  fixed  time  become  periodic  in  s. 

The  output  variables  can  then  be  Fourier  analyzed.  As  we  shall  show 
below,  the  Fourier  coefficients  represent  an  average  of  the  output  variables 
over  the  uncertainties  of  all_  the  parameters.  A  unique  correspondence 
between  the  Fourier  coefficients  for  the  frequency  and  all  Its  harmonics, 
and  the  sensitivity  of  the  output  variables  to  the  t'th  parameter  is  estab¬ 
lished.  We  compress  all  this  Information  into  partial  variances  S  which 
are  the  normalized  sums  of  the  squares  of  the  Fourier  coefficients  of  the 

fundamental  frequency  and  all  its  harmonics.  If  <  $u  for  a  given  output 

^  j 

variable,  then  this  output  variable  is  less  sensitive  to  the  t'th  parameter 
than  to  the  j'th  parameter.  Thus,  the  partial  variances  measure  the  average 
sensitivity  of  an  output  function  to  the  variation  (or  uncertainty)  of  a 
particular  parameter.  This  average  Is  over  the  range  of  uncertainties  of 
all  the  parameters,  with  their  appropriate  probability  distributions,  with 
the  exception  of  the  parameter  being  considered.  For  this  parameter,  the 
statistical  property  calculated  is  the  variance- 
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The  sensitivity  analysis  presented  here  is  nonlinear  so  that  it  permits 


us  to  examine  large  deviations  from  the  nominal  parameter  values.  In 
addition,  since  all  parameters  are  varied  simultaneously,  one  explores 
regions  of  parameter  space  where  more  than  one  parameter  is  far  from  its 
nominal  value.  Because  of  this  thorough  and  systemaci'*  exploration  of  the 
parameter  space,  it  often  turns  out  (see  e.g.  the  examples  of  Section  IV) 
that  sensitivities  of  an  unexpected  nature  are  revealed.  A  careful  study 
of  the  model  will  then  reveal  some  complex  coupling  between  variables, 
unexpected  prior  to  the  analysis,  which  leads  to  the  observed  sensitivity. 
In  this  fashion,  one  obtains  deeper  insights  into  the  structure  of  the 
complex  system  being  studied.  Another  frequent  and  important  finding  is 
that  a  number  of  sensitivity  coefficients  corresponding  to  a  large  set  of 
parameters  turn  out  to  be  negligible.  This  permits  one  to  reduce  the 
complexity  of  the  set  of  model  equations  and  focus  one's  attention  on  a 
greatly  reduced  set  of  equations.  The  applications  that  we  present  in 
Section  IV  will  exhibit  this  feature. 

The  body  of  this  paper  is  organized  Into  four  sections.  In  Section 

II  we  develop  the  theory  of  our  method  of  sensitivity  analysis.  Section 

III  details  the  implementation  of  the  method  and  in  Section  IV  we  present 
several  examples  of  its  use  to  date.  In  Section  V  we  list  a  number  of 
unsolved  problems  on  which  further  research  would  be  most  useful. 

As  previously  mentioned,  the  implementation  of  the  method  is  possible 
without  a  knowledge  of  the  details  of  the  analytic  development.  Readers 
primarily  interested  in  the  application  of  sensitivity  analysis  to  some 
specific  problems  can  by-pass  Section  II  and  go  directly  to  Section  III. 
The  examples  of  Section  IV  will  serve  to  show  the  capability  and  range  of 
applicability  of  the  method. 
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The  analytic  development  of  Section  II  is  divided  into  several  sub¬ 
sections  which  correspond  to  the  separate  steps  that  make  up  the  sensitivity 
analysis.  We  first  introduce  the  parameter  uncertainties  (or  variations) 
together  with  a  search  curve  procedure  for  exploring  the  values  of  the  out¬ 
put  function  at  different  points  in  parameter  space.  The  path  of  the  search 
curve,  parametized  by  a  search  variable  s  ,  is  closed  on  Itself  in  the 
parameter  space  so  that  the  output  function  is  periodic  in  the  search 
variable  and  can  be  Fourier  analyzed.  Since  the  search  curve  is  closed  it 
cannot  cover  all  points  In  parameter  space.  The  interpretation  of  the 
Fourier  coefficients  as  measures  of  sensitivity  is  exact  only  if  the  search 
curve  were  to  cover  all  the  parameter  space.  To  explore  this  error,  we 
Introduce  another  space,  "e-space",  which  has  the  same  dimension  as  the 
parameter  space,  and  find  that  the  error  can  be  estimated  and  controlled. 
Next  we  note  that  the  Fourier  coefficients  must  be  evaluated  numerically 
as  finite  sums  and  we  take  into  account  the  errors  introduced  by  this  pro¬ 
cedure.  The  Fourier  coefficients  are  then  combined  in  a  well  defined 
prescription  to  yield  the  desired  sensitivity  coefficients  and  the  precise 
meaning  of  these  coefficients  is  detailed.  This  Is  followed  by 
a  discussion  of  our  original  procedure^  and  the  relation  of  our  nonlinear 
sensitivity  analysis  with  the  linear  types  of  sensitivity  analysis  usually 
employed  in  the  past. 
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II.  Fourier  Method  of  Sensitivity  Analysis 

A.  Parameter  Probability  Distribution 

We  assume  that  a  mathematical  model  of  the  desired  physical  system  has 
been  constructed.  For  many  problems  the  equation  set  takes  the  form 

d  f(t) 

— 3t—  *£.(!.k)  f(t=0)  =  ^  (2.1) 

with  f  a  vector  of  m  output  functions  (f^ ,  fg,  ...  fm)  at  time  t,  ]c  a  vector 

of  n  parameters  (k^ ,  kg,  ...  kn)  which  couple  the  nonlinear  equations 

represented  by  F,  and  with  f(t*0)  =  (f^o),  fg(o) . fm(o))  a  vector 

of  given  initial  conditions.  We  assume  that  this  set  of  equations  can 

be  solved  numerically  for  f(t)  for  any  time  t  once  f_  and  k  are  specified. 

Now  consider  each  of  the  parameters  k  to  range  over  a  continuous 

4 

set  of  values  k^  <  ®.  We  write 

*  g4(ue)  -“<uA<  +  »,  4*1,2 . n  (2.2) 

where  g^u^  denotes  some  function  of  u^.  The  variable  ua 

serves  to  vary  the  kr  In  references  (1-3)  we  chose  the  particular  form 

k4  *  *4°^  exP  w1th  *4°^  the  nominal  value  of  the  parameter  to  permit  a  wide 

range  of  variation  of  k  with  u.  The  u^  are  now  considered  to  be  independent  random  variables 

with  their  respective  probability  densities  P  (u  )  such  that  P  (u  )du  gives 

the  probability  that  the  random  variable  u  lies  in  the  range  (u  ,  u  +  du  ). 

*  £  £  £ 

Since  the  variables  u£  are  not  correlated,  then  the  probability  density  P(u) 
is  given  the  product  of  the  Pfu^),  i.e. 
n 

P(u)  *  n  P  (u  )  *  (2.3) 

4*1  *  * 
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By  way  of  example,  for  many  problems  these  densities  will  be  more  or  less 

normally  distributed  with  widths  dependent  on  one's  knowledge  of  the 
dispersion  of  the  values  of  the  parameters  k^  around  the  nominal  value  k^. 

We  now  introduce  averages  over  these  densities  which,  as  we  shall  show, 
will  lead  to  useful  measures  of  sensitivity.  The  n  dimensional  u^space 
average  of  a  function  f(iO  fs  defined  as 

<f (u_)>  =  /du^  f (uj  P(uJ.  (2.4) 

The  "brute  force"  method  of  sensitivity  analysis  corresponds  to  picking 
a  grid  of  points  in  the  u_-space,  evaluating  the  solution  of  the  equation 
set  (2.1)  for  the  values  of  the  parameters  k/u)  at  these  grid  points  and 
examining  how  the  output  functions  change  with  k.(u). 
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B.  The  Search  Curve 


In  our  method  we  construct  a  search  curve  in  the  u-space  which  is  a 
path  in  u^space  parametized  by  the  search  variable  s.  The  search  curve  is 
constructed  so  that  the  n-dimensional  u^-space  average  of  the  output  function, 
Eq.  (2.4)  equals  the  one  dimensional  s-space  average. 

The  translation  of  the  probability  density  of  the  u^'s  into  s-space 
is  made  by  the  introduction  of  n  transformation  functions  ( 2.=!  ,2,. .  .n) 
such  that 

u^  =  G^sin-^Sy  (2.5) 

In  these  equations,  the  are  a  set  of  incommensurate  frequencies, 
with  one  frequency  assigned  arbitrarily  to  each  component  u^of  the  vector  u_. 
The  set  of  equations  (2.5)  then  is  a  parametric  representation  of  an  n-dimen¬ 
sional  curve  in  the  space  of  the  variables  (u-j ,  u 2,  ...  un).  Variation  of 
s  over  the  range  —  s  s  s  +  °°  generates  a  curve  which  traverses  this  n-dimen¬ 
sional  parameter  space  infinitely  often  in  each  direction,  but  with  a  relative 
rate  of  traversal  in  each  direction  which  is  proportional  to  the  frequency 
assigned  to  that  direction. 

Within  this  broad  statement,  the  detailed  shape  of  the  curve  depends 
upon  the  specific  functional  forms  chosen  for  the  transformation  functions 
G  .  Inasmuch  as  we  wish  to  obtain  a  specific  distribution  P£(u£)  in  the 
2-th  direction,  it  is  necessary  to  choose  G^  so  that  the  fraction  of  the 
total  arc  length  of  the  curve  which  lies  between  the  values  u  and  u  +  du 

X»  )b  M 

is  equal  to  P^u^Jdu  .  The  conditions  for  this  equivalence  were  deduced 
(4) 

by  Weyr  '  and  in  terms  appropriate  for  a  sensitivity  analysis,  derived  in 
an  earlier  paper  in  this  series. ^  It  was  found  there  that  if  G  (x  )  is 

Xi  JL 


taken  to  be  the  solution  of 


f 


*0  -  X2)1/2  p4(g£) 


dG.(x) 


dx 


=  1 


(2.6) 


with  the  boundary  condition  Gt(0)  «  0,  then  the  arc  length  along  this  curve 
is  distributed  in  accordance  with  the  previously  stated  requirement.  This 
equation  can  be  solved  by  quadrature  for  any.  distribution  function  P^.  Thus, 
it  is  possible  to  construct  a  one-dimensional  manifold  which  covers  para¬ 
meter  space  exactly  in  accordance  with  the  requirements  of  the  statistics. 

The  curve  generated  by  the  set  of  equations  (2.5)  and  (2.6)  can  be 
referred  to  as  a  "search  curve".  As  demonstrated  by  Weyl ,  the  fact  that 
the  frequencies  u>4  are  incommensurate  guarantees  that  the  curve  is  space¬ 
filling,  by  which  we  mean  that  it  passes  arbitrarily  close  to  any  preselected 

n 

point  in  parameter  space  for  which  the  joint  distribution  n  ^(u^)^ 
does  not  vanish.  Thus,  the  search  curve  is  an  ideal  sampler  of  parameter 
space,  since  it  seeks  out,  in  the  sense  just  defined,  each  point  in  the 
space,  and  it  fills  the  space  with  a  relative  density  matched  to  the  joint 
distribution  of  the  parameters. 

The  search  curve  just  described  is  an  ideal  which  cannot  be  numer-  -.ally 
realized.  It  is  clear  that  the  frequencies  which  are  used  in  the  com¬ 
putational  analysis  cannot  be  incommensurate.  Thus,  we  are  limited  to 
commensurate  frequencies  which  must  be  chosen  to  mimic  as  closely  as  possible 
the  above  requirement  of  a  space-filling  curve.  We  shall,  in  the  sections 

that  follow,  discuss  the  errors  Introduced  by  the  use  of  commensurate  fre- 

fl  2  3l 

quencles.  As  in  our  earlier  work'  *  ’  '  we  develop  the  theory  here  for 
Integer  frequencies.  The  use  of  such  frequencies  leads  to  a  search  curve 
that  cannot  fill  the  space  but  yields  a  closed  path  in  the  u.-space.  The 
entire  traversal  of  the  u-space  Is  then  accomplished  by  letting  the  search 
parameter  s  range  between  0  and  2n. 
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C.  Fourier  Coefficients 


The  use  of  integer  frequencies  in  the  transformation  functions, 

Eq.  (2.5),  implies  that  the  u£'s  are  periodic  in  s  on  the  interval  (0,  2*). 
Since  the  output  functions  f.  are  a  function  of  s  through  the  u^'s  they  are, 
as  a  function  of  s,  periodic  on  (0,  2rr),  i.e. 


ff  (s)  =  f|  (s  +  2tt). 


(2.7) 


The  output  functions  can  thus  be  Fourier  analyzed  to  obtain  their  Fourier 
coefficients 


A 


0) .  L  r’ 

p“l  "  'o 


cos(pUjs)  f. (s)  ds; 


P  -  0,  1,  ... 


.0)  -1  f  sin(pw  s)  f. (s)  ds;  .  p  *  1 ,  2,  ... 

pUji  ir  0 

In  terms  of  an  average  over  s-space,  the  Eqs.  (2.8)  can  be  re-written  as 


=  <cos(pu.As)  fi(s)> 

(2.9) 

?Bp^  3  <sin(PV}  fi(s)> 

The  importance  of  the  Fourier  coefficients  in  Eqs.  (2.8)  is  that  if 
and  are  zero  for  all  p=l,2...  then  the  i'th  output  f4  is  Insensi- 
tive  to  the  values  of  the  Ji'th  parameter  k^.  The  Fourier  coef-  ^ 

ficients  of  the  £ 1 th  fundamental  and  all  harmonics  of  this  i'th  fundamental 
measure  the  sensitivity  of  the  output  to  the  &'th  parameter.  These  statements 
will  be  justified  in  the  following  section. 

The  Fourier  coefficients  can  be  interpreted  as  sensitivity  measures 
if  we  can  show  that  the  Fourier  coefficients  of  a  given  fundamental  and  all 
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f 


Its  harmonics  segregate  the  effects  of  each  parameter  uncertainty  on  the 

output  functions.  That  is,  if  the  Fourier  coefficients  A_  and  EL  (p=l,2...) 

Pw£ 

are  affected  by  the  uncertainty,  i.e.  range  of  variation,  in  the  i'th 
parameter  k  and  are  not  affected  by  uncertainties  in  any  of  the  other 
parameters,  then  these  Fourier  coefficients  isolate,  one  by  one,  the  un¬ 
certainties  of  the  parameters  k^  on  the  outputs.  For  the  integer  frequencies 
that  we  use  this  is  not  strictly  the  case,  but  it  is  approximately  true. 

It  is  only  for  incommensurate  frequencies  i.e.,  those  for  which 


n 

L  r.u).  f  0  -»  <  r4  <  +  oo  (2.10) 

1=1  1  1  1 


with  r..  integer,  that  the  segregation  of  sensitivity  is  exact.  However, 
we  shall  show  that  the  error  made  in  the  use  of  integer  frequencies  can  be 
quantitatively  predicted  and  controlled.  We  are  therefore  still  able  to 
use  this  very  useful  interpretation  of  the  Fourier  coefficients  as  proper 
measures  of  sensitivity. 


0.  The  e_-Space 


In  order  to  demonstrate  the  validity  of  this  interpretation  of  the 
Fourier  coefficients,  it  is  necessary  to  visualize  the  search  curve  in  an 
n-dimensional  space  which  is  different  from  the  u-space.  Let  us  introduce 
an  n-dimensional  ^-space  with  orthogonal  axes  0-|,  @2.  •••  en  defined  by 

*  oj^s  (mod  2ir ) ;  1*1,2,  ...  n.  (2.11) 

In  this  n-dimensional  space  the  search  curve  consists  of  a  series  of  parallel 
lines  with  the  separation  between  the  lines  determined  by  the  choice  of 
integers  As  simple  examples,  consider  the  two  different  choices  of 
integer  frequencies 

<4>1  =  1 ,  Wg  =  2  (2.12a) 

«1  *  11,  *z  »  13.  (2.12b) 

The  two  dimensional  ^-space  and  the  search  curves  for  these  frequency 
choices  are  shown  in  Figures  2  and  3. 

It  is  intuitively  clear  that  as  the  search  curve  does  a  better  and 
better  job  of  covering  the  ^-space,  an  integral  of  a  function  over  s-space 
can  equally  well  be  carried  out  as  a  multidimensional  integral  over  ^-space. 
The  way  to  obtain  a  uniformly  dense  coverage  of  the  0-space  is  to  choose 
w.|  and  to  be  incommensurate,  that  is  choose  and  such  that 

rl"l  +  r2w2  *  0  --<r1,r2<  +  “  (2.13) 

(4) 

for  any  choices  of  the  integers  r1  and  r2-  Integer  frequencies  cannot 
be  incommensurate  but  by  choosing  w-j  and  ^  appropriately  the  values 
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of  rj  and  which  lead  to  equality  In  Eq.  (2.13)  become  very  large.  In 
the  two  examples  we  have  chosen,  *  2  and  rg  *  -1  lead  to  equality  In 
Eq.  (2.13)  for  the  first  choice  of  Integer  frequencies  while  r-j  *  13  and 
Tg  *  -11  are  required  for  the  second  choice  of  frequencies.  Thus  the  latter 
frequency  choice  leads  to  a  better  coverage  of  the  e-space  as  is  evident 
from  Figures  2  and  3. 

So  far  we  have  focussed  on  the  search  curve's  coverage  of  e^-space  without 
reference  to  the  function  which  is  being  integrated.  The  accuracy  of 
replacing  the  s-space  average  by  a  e-space  average  not  only  depends  on  this 
coverage  but  also  on  the  values  of  the  integrand  since  we  want  to  equate 
.ff(s)P(s)ds  with  J*f(e)P(e)de_.  As  an  extreme  example,  if  the  function  f(s) 
were  constant  as  s  ranged  over  (0,  2 it)  ,  then  one  could  equate  the  s  and 
£-space  average  without  error  for  any  choice  of  search  curve.  If  f(s)  is 
a  slowly  varying  function  of  s,  then  it  will  also  be  a  slowly  varying 
function  of  £  in  £-space  and  the  error  in  equating  the  s  and  £-space 
integrations  will  be  small  for  any  choice  of  search  curve.  This  error  is 
likely  to  increase  when  f(s)  and  f(e)  are  rapidly  varying  functions  over 
the  ranges  of  s  and  a.  In  equating  s  and  e-space  averages  one  must  therefore 
consider  both  the  coverage  of  space  by  the  search  curve  and  the  "smoothness" 
property  of  the  output  function  f(s). 


E.  Relation  Between  s  and  a-Space  Averages 

Since  the  interpretation  of  the  Fourier  amplitudes  as  sensitivity 
coefficients  will  be  made  in  e_-space  it  is  necessary  to  relate  quanti¬ 
tatively  the  s  and  e_  space  averages.  This  relation  is  obtained  by 
expressing  the  output  function  f(s)  as  a  function  of  e.  and  noting  that,  by 
construction,  f(ej  is  multiply  periodic  in  e_  on  (0,  2v).  We  may  thus  expand 
f(ej  in  a  multiple  Fourier  series 

f(i>  -  I  l  l  cr1r,...r„ex',[-((slrlte2r2t--'f9nrr>] 

rl  r2  rn  1  2  " 

*  l  cr  exp[-ie/r]  ,  -  »  <  r^  <  +®  (2.14) 

The  s-space  Fourier  coefficients  defined  in  Eqs.  (2.9)  are  a  subset 
of  all  possible  such  coefficients,  namely  the  subset  corresponding  to  a 
fundamental  frequency  u  and  all  its  harmonics  pw4(ps2,3,...).  In  e^- space 
we  need  the  corresponding  Fourier  coefficients  for  our  sensitivity  analysis. 
As  discussed  previously,  these  s  and  £-space  coefficients  are  not  equal 
because  the  search  curve  does  not,  for  a  closed  path,  cover  the  entire 
e-space.  The  difference  between  these  Fourier  amplitudes,  i.e. 

27  f  f(s)e1pVds  -(27)  /•••/def(i)exp[ipeJl] 

is  the  difference  between  an  Integral  calculated  with  the  search  curve,  a 
line  through  ^  space,  and  an  integral  calculated  over  the  entire  e_  space. 
Written  as  averages,  this  difference  is 
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(2.15) 


<f(s)eifV>  -  <f(e>ipet>  =  e<f(0_)eip04> 

where  e  measures  the  relative  error  In  equating  the  s  and  0_-space  Integrals. 
Equivalently,  we  can  write 


Cpui  '  c000...pr..000  *  ec000...p  ... 000 

where  C.  is  the  complex  Fourier  coefficient  defined  by 

J 


(2.16) 


cj  *  r  (aj  -  1Bj» 

(2.17) 

c-j  "  cj*  *  7  +  iBj> 

and  where  Cqq  p  0Q  is  defined  in  Eq.  (2.14). 

In  reference  (3)  and  in  Appendix  I  of  this  paper  we  show  that  by 
appropriate  choices  of  frequencies  u  ,  this  error  can  be  made  as  small  as 

X* 

desired.  Qualitatively,  as  the  frequencies  become  more  incommensurate  the 
coverage  of  the  ^-space  (and  u-space)  by  the  search  curve  improves  and  the 
error  made  in  equating  s  and  e.-space  averages  decreases.  This  can  be 
quantified  through  the  introduction  of  an  index  M  which  is  the  order  of 
interference.  Interferences  arise  when  frequencies  in  the  original  set  ^ 
can  combine  to  form  another  frequency  of  the  set.  Thus,  if  we  have  three 
parameters  with  associated  frequencies  ^ ^  and  u>3  and  if,  for  instance, 
w-j  +  uig  *  W3,  an  interference  has  occurred.  As  M  increases  the  interferences 
are  postponed  to  higher  and  higher  harmonics  and  the  error  decreases. 

Assuming  that  e  is  negligible  by  appropriate  choice  of  the  set  ^  we 
then  obtain 

Cpu>.  "  c000...p  ...000  "  cp,  t2,18) 

1  1  1 
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as  the  desired  equality  of  s  and  e^-space  averages.  The  Fourier  coefficient 
Cqqq  p  qqq  reflects  the  sensitivity  of  the  output  function  to  the 
uncertainty  of  the  t,'th  parameter  (and  only  to  the  Jt'th  parameter)  since 


all 


r4  *  0  for  all  i  f  a  in  the  e-space  Fourier -coefficient  c  .If 

i  -  rl  r2  “  *  rn 

the  Fourier  coefficients  Cqq  q00  are  zero  for  p4  *  1,2*...,  then 

we  can  conclude  that  the  output  function  is  not  sensitive  to  the  parameter  k 

Since,  by  (2.18),  the  e-space  Fourier  coefficients  c_  equal  the  s-space 

-  Pi 

coefficients  C„  ,  we  conclude  that  C„  is  an  appropriate  measure  of  the 

P“i 

sensitivity  of  the  output  to  the  i'th  parameter. 

The  conclusion  that  the  Fourier  coefficients  C  measure  just  the 
sensitivity  to  k  ,  can  also  be  verified  via  an  s-space  analysis.  The  output 

jL 

as  a  function  of  s,  f(s),  is  constructed  by  expressing  each  parameter  u^ 
as  a  periodic  function  of  s  with  frequency  <*>.  (i=l,2,...n)  as  in  Eq.  (2.5). 
Thus,  f(s)  consists  of  a  sum  of  terms  which  oscillate  at  all  possible  linear 
combinations  of  the  n  frequencies  a.  If  we  single  out  the  frequency  pu 
via  the  Fourier  coefficient  Cpw  ,  and  if  the  frequencies  are  incommensurate, 
no  linear  combination  of  the  other  frequencies  can  add  up  to  form  and 
interfere  with  the  effect  on  the  output  from  the  n'th  parameter  k  .  Since 
we  do  use  integer  frequencies,  it  is  necessary  to  introduce  the  concept  of 
order  of  interference  M  as  defined  in  Eq.  (A1.4).  Thus,  as  long  as  we 
restrict  our  attention  to  Fourier  coefficients  C„  for  which  there  are  no 
interferences,  or  for  which  the  interferences  have  been  "postponed",  the 
sensitivity  is  due  entirely  to  the  uncertainty  In  k  . 
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F.  Finite  Fourier  Coefficients 

In  performing  the  Fourier  analysis  on  the  computer,  the  s-space  inte¬ 
gration  must  be  approximated  by  a  finite  summation.  This  procedure  introduces 
a  further  error  into  the  method,  which  must' be  analyzed.  It  is  important 
to  minimize  this  error  as  efficiently  as  that  due  to  interferences.  The  finite 
sums  are  obtained  by  taking  points  from  the  search  curve  in  e^space.  If 
we  use  a  large  number  N  of  points  for  the  summation  we  will  obtain  an 
excellent  approximation  to  the  s-space  integral.  However,  if  the  order  of 
interference  M  is  low,  the  accuracy  of  the  interpretation  of  the  s-space 
average  as  a  sensitivity  measure  Is  not  good.  For  example,  if  we  fix  N  at 
100,  in  the  2-dimensional  examples  discussed  in  section  C  above,  the  first 
search  curve  leads  to  the  set  of  points  shown  in  Figure  2  while 
the  Second  search  curve  leads  to  the  set  of  points  shown  In 
Figure  3.  Even  though  the  first  search  curve  is  better  approximated  by  the 
100  points  quadrature  than  the  second  one  (since  the  density  of  points  Is 
higher  for  the  first  curve),  the  overall  distribution  of  points  In  the 
second  case  is  far  more  uniform.  We  therefore  would  obtain  a  more  accurate 
sensitivity  measure  in  the  second  case  than  in  the  first. 

The  error  introduced  via  the  use  of  finite  sums  for  the  s-space  Inte¬ 
grations  can  be  analyzed  in  much  the  same  fashion  as  the  interference  problem. 
We  denote  the  numerical  approximation  to  the  exact  Fourier  coefficient 
C„  by  C*  where 

f(sq)exP[ip“tsq]  (2J9) 

with 

sq  «  ^  q  *  1 ,2, . . ,N.  (2.20) 
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Here  N  is  the  number  of  points  used  in  the  quadrature  formula  and 

(q  =  1,2,...N)  labels  the  individual  points  which  are  chosen  to  be  equally 

spaced  along  the  search  curve.*  Following  the  procedure  used  in  Eq.  (2.15), 

* 


we  now  form  the  difference  between  the  s-space  quadrature 
9j-space  average  C  to  obtain 


i 


and  the 


(2-21) 


where,  see  Eq.  (2.18),  c  s  cQ0  >i000. 

I 

The  coefficient  e*  in  Eq.  (2.21)  is  a  measure  of  the  error  introduced 
by  approximating  the  e_-space  integration  as  a  finite  sum  at  selected 
points  in  s-space.  These  errors  are  in  addition  to  those  arising  from 
interferences.  The  evaluation  of  e*.  as  a  function  of  w,  N  and  the  output 
function,  is  presented  in  Appendix  II. 

It  is  important  to  have  a  ready  estimate  of  the  error  term  when  doing 
sensitivity  analysis  without  becoming  involved  in  long  calculations  for 
each  particular  case  under  investigation.  In  response  to  this  need  we  have 
in  ref.  (3)  developed  bounds  on  the  error  which  depend  on  w,  N,  and  the  type 
of  output  functions  being  investigated.  We  present  the  salient  ideas  of 
this  method  in  Appendix  III. 


*  The  equal  spacing  is  not  necessary;  however,  the  analysis  and  control 

of  the  error  is  facilitated  by  this  choice.  The  use  of  other  spacings 

may  actually  lessen  the  error.  See  e.g.  V.I.  Krylov  and  L.G.  Kruglikova, 

Handbook  of  Numerical  Harmonic  Analysis,  Israel  Program  for  Scientific 
Translations,  LTD.  Jerusalem,  1969. 
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G.  Partial  Variances 


The  Fourier  coefficients  with  which  we  have  been  concerned  so  far 
correspond  to  the  fundamental  and  harmonics  of  each  frequency  They 
measure  the  sensitivity  of  an  output  to  the' variations  in  all  the  parameters 
such  that  if  there  is  a  sensitivity  to  the  2,'th  parameter  it  will  show  up 
only  in  the  Fourier  coefficients  C  (p=l,2,...).  The  uncertainty  in  the 
other  parameters  is  accounted  for  by  averaging  over  their  respective 
distributions.  It  is  apparent  that  other  Fourier  coefficients  contain 
additional  information  such  as  for  instance  the  joint  sensitivity  to 
parameters  k  and  k..  From  this  point  of  view,  we  can  consider  the  variation 

£  J 

of  the  output  arising  from  the  uncertainties  in  all  the  parameters  k^ , 

1=1,  2,  ....  n  and  their  couplings  to  be  characterized  by  the  variance  of 
the  output  function 

o2  =  <f2>  -  <f>2  .  (2.22) 

The  interpretation  of  this  variance  can  be  explored  by  expressing  it  in 
terms  of  the  s-space  Fourier  coefficients  Cj  (and  the  cos  and  sin  coefficients 
A.  and  B.)  by  using  the  Fourier  series  expansion  of  f(s)  in  Eq.  (2.22). 

J  J 

This  yields 

o2  -  1'  c2  *  l  l  (A2  +  B2)  .  (2.23) 

j=~  J  c  j=l  J  J 

where  the  prime  on  the  sum  excludes  the  j=0  term.  We  recognize  Eq.  (2.23) 

as  Parseval's  Theorem  of  Fourier  analysis.  However,  we  are  calculating 

*  *  * 

finite  sums  for  the  numerical  Fourier  coefficients  C ^ (Aj ,Bj )  and  must  there- 

2* 

fore  modify  Parseval's  Theorem  to  obtain  the  numerical  variance  o  .  This  is 
just 
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2* 

The  variance  0  is  thus  seen  to  be  the  sum  of  the  squares  of  the  Fourier 
coefficients  of  all  integer  frequencies  which  enter  into  an  N  point 
quadrature  formula. 

2 

We  now  construct  the  part  of  o  that  corresponds  to  the  variance  of 

the  output  arising  from  the  i'th  parameter  uncertainty  when  the  output  is 

averaged  over  the  uncertainties  in  all  other  parameters.  To  do  this,  we 

first  integrate  f(k)  and  f  (k_)  over  all  uncertainties  in  the  parameters 

with  the  exception  of  the  i'th  parameter.  We  then  calculate  the  variance 

ot  for  these  partially  integrated  output  functions  using  equations  (2.22) 

and  (2.23).  The  details  of  this  calculation  will  be  found  in  Appendix  IV. 

2* 

This  result  is  then  modified  to  yield  the  numerical  variance  corresponding 

2  * 

to  the  finite  sums  of  Eq.  (2.24).  The  ratio  of  this  variance  to  the 

total  variance  o  of  Eq.  (2.24)  is  denoted  by  S  and  is  the  partial 

I 

variance,  i.e. 


NZ2  |  * 
p«-(N/2-iyCpu)* 


T 

j=- (N/2-1 ) 


★ 

co 


(2.25) 


* 

The  partial  variance  S  is  the  fraction  of  the  variance  of  the  output  function 
due  to  the  uncertainty  in  the  parameter  k4  when  the  output  function  is 
averaged  over  the  uncertainties,  and  coupling  of  uncertainties  of  all  the 
other  parameters  k^ ,  i£t. 

The  partial  variance  can  thus  serve  as  an  excellent  measure  of  the 


sensitivity  of  the  output  to  the  uncertainty  of  the  t'th  parameter.  It  can 
be  used  to  compare  quantitatively  sensitivities  due  to  uncertainties  in  the 


different  parameters  k^ ,  i=l,2,...n,  since  all  variances  are  computed 

relative  to  the  total  variance  o  .  Clearly,  the  smaller  the  partial 

variance,  the  less  the  effect  of  the  changes  of  k  on  the  output  function. 

★ 

We  can  therefore  order  the  S  's  as  a  function  of  i  to  obtain  an  ordered 

“2 

list  of  sensitivities  for  a  given  output  function  and  between  different 
output  functions. 

★ 

The  partial  variances  S  will  have  the  above  interpretation  only  if 

“2  (see  Appendix  II) 

the  interferences  and  aliasing  difficulties^are  minimized  by  proper  choice 

★ 

of  at,  and  N.  Since  S  involves  Fourier  coefficients  of  high  index,  the 
appropriate  values  of  w,  and  N  for  the  minimization  of  these  errors  unfor¬ 
tunately  involve  prohibitively  large  numbers  of  function  evaluations. 

However,  as  discussed  previously,  the  Fourier  coefficients  decrease  in 
amplitude  as  their  index  increases.  As  a  practical,  computational  matter, 
our  experience  has  shown  that  only  a  few  harmonics  of  a  fundamental  need 

be  calculated  before  the  amplitude  becomes  negligible.  Thus,  the  relative 

★ 

errors  in  comparing  the  partial  variances  S  for  different  parameters  k. 

W2  1 

are  small  even  when  only  the  Fourier  amplitudes  of  the  fundamental  and  the 
first  few  harmonics  are  utilized. 

Let  us  return  briefly  to  the  unstarred,  i.e.  analytic  values  of  the 
2 

variance  o  (Eq.  2.23)  and  the  corresponding  partial  variances  S  given  by 

°°  2  1 
2  Tcn 

S  ■  (2.26) 

2  0  .  z  c  i 

2 

where  the  prime  on  the  sums  excludes  the  p=0,  j=0  terms.  Note  that  ol  in 

the  numerator  involves  only  the  sum  of  the  squares  of  the  Fourier  coefficients 

of  the  fundamental  and  all  harmonics  of  the  2'th  frequency  while  the 
2 

variances  o  in  the  denominator  is  constructed  as  the  sum  of  squares  of  the 


Fourier  coefficients  of  al 1  the  integer  frequencies.  It  is  readily  apparent 

n 

from  Eq.  (2.261  and  the  development  below  why  the  sum  7  S  of  the  partial 

t=l  ui 

variances  will  not  equal  unity.  Writing  the  variance  o2  in  terms  of  the 
c_-space  Fourier  coefficients  according  to  Parseval's  theorem  yields 


2 

0 


cd  oc  or 


■  v  v  ...  ric(Pl,p2 . pji 

Pl  p2  Pn 


2 


(2.27) 


where  the  prime  indicates  that  the  term  with  all  p^'s  equal  to  zero  is 
omitted.  It  is  suggestive  to  rearrange  the  sum  in  Eq.  (2.27)  into  groups 
of  terms  where  successively  larger  subgroups  of  the  coefficients  Pj,P2,...,p 
are  nonzero.  We  define 


9  9 

o  =  l  jc(o,...,p  ....  ,0 )  I 
x  p.=l  1 

2CC  00  _ 

—  9 

o..  *  _  £  Jc(o,...,p,...,p.,...,o)j 

J  pt.l  Pj.l  J 


2  2 

=  £  I  I  I  c(o , .  .  .  ,p.  , .  .  .  ,p  • , .  .  .  ,P|^  , .  .  .0 )  | 

4Jk  P£=l  p.=l  pk-l  *  3  k 

etc. 


so  that 


(2.28a) 

(2.28b) 

(2.28c) 


,  n  ,  n  n-1  ?  n  n-1  n-2  9 

o2  ■  I  -  ♦  I  I  As*  l  l  l  '  Jk  * 

£=1  *■  i  =  l  j  =  l  8  =  1  i  =  l  k  =  l  X'JK 


1-1  j=l  k=l 


(2.29) 


is  written  as  a  sum  of  terms  with  successively  more  complex  contributions 
to  the  total  variance  o  . 

We  have  shown  above  that  the  first  term  of  the  decomposition  (2.29) 
corresponds  to  the  part  of  total  variance  arising  from  the  t'th  parameter 
uncertainty  when  the  output  is  averaged  over  the  uncertainties  in  all  other 
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rate  coefficients.  By  the  same  methods  one  can  also  construct  by 

JCJ 

integrating  f(k)  over  all  but  the  two  parameters  k^  and  kj  to  obtain  ^(^.kj) 

2 

and  then  forming  the  variance  o  .  of  this  function  via  Eq.  (2.22).  One 

*  J 

can  proceed  in  an  analogous  fashion  to  generate  the  other,  more  complex, 
variances  in  Eq.  (2.29).  These  higher  partial  variances  are  thus  seen  to 
contain  increasingly  more  detailed  information  about  the  coupling  of 
sensitivity  among  larger  and  larger  groups  of  parameter  uncertainties. 

In  our  own  investigations  to  date  we  have  not  exploited  these  higher 
partial  variances  to  obtain  more  detailed  information  about  the  sensitivity 
of  our  test  systems  but  we  hope  to  explore  their  properties  in  some  subsequent 
studies. 
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H.  Fundamental  Fourier  Coefficients 

In  the  original  formulation  of  sensitivity  analysis,  references  (1) 
through  (3),  our  efforts  centered  on  the  Fourier  sine  coefficients  of  the 
fundamental  frequencies,  the  in  our  present  notation.  As  shown  there, 
these  Fourier  coefficients  can  be  expressed  as  £-space  averages: 


B 

“l  a*  X 


(2.30) 


where  the  average  was  taken  over  the  u-space  probability  density 


n 

P(u)  =  n  a  ./cosh{a  .u.)  (2.31) 

j=l  3  J  J 

with  the  P(u.,a.)  being  bell  shaped  distributions  whose  widths  are  deter- 

J  J 

mined  by  the  parameter  aj.  We  have  in  these  papers  also  explored  various 
other  probability  densities. 

The  Fourier  coefficients  Bw  are  thus  seen  to  be  directly  related  to 
the  rate  of  change  of  the  output  function  with  respect  to  the  *'th  parameter, 
averaged  over  the  uncertainty  in  all  the  parameters.  This  appealing 
interpretation  of  the  fundamentals  does  not  yield  as  sharp  a  sensitivity 
measure  as  do  the  partial  variances  since  the  integrand  of  the  average 
<df/du£>  is  not  necessarily  a  positive  function  and  consequently  could 
erroneously  indicate  a  small  sensitivity  by  fortuitous  cancellation  in 
different  regions  of  parameter  space.  The  partial  variances  are  not  subject 
to  this  difficulty;  if  S,  is  zero  or  smaller  it  is  definitely  due  to  lack 

i 

of  sensitivity  of  the  output  to  the  *'th  parameter. 
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I.  Relation  to  Linear  Sensitivity  Analysis 

We  discuss  here  the  relationship  of  our  nonlinear  method  of  sensitivity 
analysis  to  the  more  usual  methods  of  linear  sensitivity  analysis.  In  the 
linear  methods  one  computes,  one  way  or  another,  the  derivatives  dc/du^^g 

for  t=l,2 . n.  This  validity  of  this  procedure  must  certainly  be  suspect 

for  output  functions  which  deviate  significantly  from  linearity  in  the 
uncertainties  in  one  or  more  of  the  parameters.  It  is  thus  clearly  of  interest 
to  investigate  the  conditions  under  which  our  Fourier  amplitude  analysis 
reduces  to  the  linear  analysis. 

To  demonstrate  this  relationship,  we  express  the  average 

<f(x)>  h  f  •••  /  dx  p(x)f(x)  (2.32) 

of  a  function  f(x)  in  terms  of  the  Taylor  series  expansion  of  f(xj.  The 
multivariate  probability  distribution  function  p(x)  can  be  written  formally 
as 


p(x)  = 


m  r  T 

lli} 


(2.33) 


/  u .  \ 

where  6V  1 y (x^ )  is  the  w^'th  derivative  of  the  Dirac  delta  function  with 
respect  to  its  argument. 


OD  oo 


and  where 


(2.34) 


are  the  multivariate  moments  of  p(x).  Using  this  expression  for  p(x)  in 
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Eq.  (2.32)  yields 


(2.35) 
x=0 

This  is  the  desired  relationship  between  the  average  of  a  function  of  a  set 
of  variables  and  the  coefficients  of  the  Taylor  series  expansion  of  the 
function. 

This  analysis  is  readily  applicable  to  the  probability  distributions 
and  output  functions  discussed  in  the  previous  sections.  Let  us  choose, 
for  instance,  the  probability  density 

n 

P(u)  =  n  a./cosh  (a.u.) 
j*l  3  33 

of  eq.  (2.31)  with  moments 


where  is  the  Euler  number  of  index  K.  The  output  function  In  u^space 
which  we  will  consider  is 

f<H)  *  r  Tu**-  (2-37) 

£  l 

since.by  Eq.  (2.30),  the  average  of  this  function  with  the  probability  density 
P(u)  of  Eq.  (2.31)  yields  the  Fourier  coefficient  B  ,  i.e. 


-28- 


where  the  prime  on  the  y  sum  excludes  the  ^  =  ^  =  • • •  *  un  =  0  term  . 

In  order  for  linear  sensitivity  analysis  to  be  valid,  the  higher 
order  terms  of  Eq.  (2.4Q)  must  be  small  compared  to  the  first,  i.e.  the 
linear,  term.  These  higher  order  terms  will  be  small  only  if  the  output 
function  is  essentially  linear  in  the  uncertainties  of  all  the  parameters. 

The  validity  of  a  linear  sensitivity  analysis  can  therefore  not  be  established 
until  some  evaluation  or  estimation  of  these  higher  order  terms  has  been 
carried  out.  Our  Fourier  amplitudes,  as  shown  in  Eq.  (2.40),  do  include  the 
effects  of  the  higher  order  terms  and  therefore  represent  nonlinear  sensitivity 
measures. 

The  results  of  this  section  can  easily  be  transcribed  to  show  that  the 

partial  variances,  S  ,  as  defined  in  Eqs.  (2.25)  and  (2.26),  are  also  non- 

i 

linear  sensitivity  measures  which  include  the  effects  of  the  higher  terms  of 
the  expansion  (2.40).  This  Is  readily  evident  from  the  definition  of  the 
variances  o  and  in  terms  of  the  Fourier  amplitudes,  Eqs.  (2.23)  and  (A4.2). 
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III.  Implementation 

The  sensitivity  coefficients  are  given  by  sums  of  the  solutions  of 
the  equation  set  at  selected  points  in  the  uncertainty  space.  From  the 
point  of  view  of  computation,  all  that  is  required  is  the  solutions  of  the 
equation  set  for  different  combinations  of  the  parameters  (k^  ,kg » - . . kn ) . 

The  analysis  of  section  II  shows  which  combinations  of  solutions  to  sum 
and  interprets  the  sum  as  a  sensitivity  coefficient. 

The  parameters  are  varied  simultaneously  by  relating  them  to  a  search 
variable  s  and  a  frequency  set  w  =  Uj  ,o>2*  •  •“»„)  as  given  in  Eqs.  (2.2)  and 
(2.5),  i.e. 


kz  *  9JtCG£(sina)jls)] 


U*1  ,2,...n). 


(3.1) 


The  frequency  «  and  transformation  functions  G  are  chosen  such  that  the 
* 

path  in  k-space  induced  by  varying  s  traverses  space  in  accordance  with 
the  probability  density  P  (u  )  representing  the  uncertainty  in  each  parameter 

X  X 

k^.  As  long  as  k^ and  s  space  averages  are  properly  connected,  sensitivity 
coefficients  can  be  calculated  as  s-space,  i.e.  one  dimensional,  quadratures. 

The  appropriate  s-space  averages  to  calculate  are  the  s-space  Fourier 
coefficients  defined  in  Eq.  (2.8).  The  analysis  of  section  II  demonstrates 
that  the  combination  of  Fourier  coefficients  which  we  call  the  partial 
variances 


I  (Ia'I1!2  *  |b<!>|2) 


SV  ‘  J  (|A]’'|2t  |Bj(1)|2) 


(3.2) 


measure  the  sensitivity  of  the  1th  variable  to  the  uncertainty  in  the  t’th 


parameter.  The  sensitivity  coefficient  '  is  that  part  of  the  variance 


of  the  ith  output  which  arises  from  the  uncertainty  in  the  &'th  parameter 
when  the  output  is  averaged  over  the  uncertainties  In  all  other  parameters 
(cf.  section  I1G). 

The  Fourier  coefficients  can  be  calculated  if  solutions  of  the  equation 
set  eq.  (2.1)  are  known  for  values  of  prescribed  by  the  transformation 
functions  6^  and  frequencies  w.  Since  we  assume  that  the  reader  has  available 
a  method  of  solution  of  the  equation  set,  all  that  need  be  done  to  obtain  the 
sensitivity  coefficient  is  to  "add11  the  solutions  according  to  Eq.  (3.2). 

Implementation  of  the  above  scheme  on  a  computer  requires  certain 
compromises  which  lead  to  approximations  to  the  sensitivity  coefficients. 

We  have  shown  in  section  II  that  these  approximations  are  controllable; 
here  we  investigate  their  impact  on  the  numerical  calculation  of  the 
sensitivity  coefficients. 

In  brief,  we  are  replacing  a  multidimensional  integral  over  the  uncer¬ 
tainty  space  (the  k_  space)  by  a  quadrature  formula  which  is  a  sum  over  N 
points  in  k  space.  The  quadrature  requirements  fall  loosely  Into  three  areas: 
1)  use  of  integer  frequencies  u;  2)  use  of  a  finite  number  of  points  N; 

3)  the  choice  of  frequencies  u.  We  now  discuss  these  issues  in  turn. 

1 .  Integer  Frequencies 

If  the  search  curve  is  to  come  arbitrarily  close  to  every  point  in 
space,  the  frequencies  that  define  this  curve  must  be  Inconmensurate: 
n 

I  r.n,  f  0  (r. 's  integer)  (3.3) 

1«1  1  1  1 

A  consequence  of  this  condition  is  that  at  most  one  of  the  frequencies 
can  be  rational,  with  all  others  being  Irrational.  But  a  computer  can  only 
represent  an  irrational  number  approximately.  Of  course,  the  difference 
between  the  Irrational  and  its  rational  approximation  can  be  made  smaller 
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and  smaller  by  resorting  to  representations  of  successively  higher  precision. 
But  eventually,  one  will  run  out  of  computer  memory,  so  that  there  is  a 
real  limit  to  the  accuracy  of  the  representation. 


Once  we  accept  the  limitations  of  rational  approximation  of  irrational 

numbers,  we  realize  that  the  frequency  set  as  represented  on  the  computer 

cannot  be  truly  irrational.*  It  therefore  becomes  necessary  to  define  a 

sequence  of  approximations  to  incommensurability.  This  was  done  in  the 
/ 1  ?3l 

previous  papers'  *  *  '  as  follows: 

A  set  of  rational  numbers  (i  =  l,2,...n)  is  approximately  incommen¬ 
surate  to  order  M  if 

1  t  0 

i=l  1  1 

for  (3.4) 

1  hi  sM  + 1 

1=1 

with  M  an  integer  at  our  disposal.  It  should  be  clear  from  this  definition 
that  incoimtensurability  corresponds  to  M  ■>». 

Henceforth  we  assume  that  the  frequencies  are  approximately  incommensurate 
to  order  M. 

If  we  are  now  dealing  with  rational  numbers  we  may  as  well  take  them 
to  be  integers.  The  correspondence  between  the  rational s  and  integers 

is  simple  to  establish  with  the  introduction  of  the  least  common 
integer  multiple:  is  defined  as  the  smallest  integer  such  that 

(z  *  l,2,...n)  (3.5) 

are  integers  for  all  z. 

The  search  curve  with  u  Integer,  Eq.  (3.1),  is  a  closed  curve  in  k^ 
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space  since  for  s  outside  the  range  (-??,*},  u^  must  repeat  a  value  from  the 
range  (-tt.tt). 

The  total  arc  length  of  the  closed  search  curve  defined  above  will 
increase  with  increasing  M,  the  order  of  incommensurability.  Thus,  eq. 

(3. T)  serves  to  define  a  sequence  of  approximate  search  curves  which 
success ively  become  more  ne'arly  space  filling  as  M  -+■  «. 

We  must  balance  the  increased  accuracy  which  one  obtains  by  using  a 
longer  search  curve  against  the  increasing  computation  time  required  to 
evaluate  the  output  functions  for  more  of  the  k_  space. 

2.  Discrete  Sampling 

The  search  curve  given  in  eq.  (3.1)  does  not  yet  completely  fix  a 
sample  of  parameter  space  which  can  be  utilized  in  real  problems  because  the 
number  of  points  on  the  search  curve  is  uncountably  infinite.  We  must 
select  a  finite  subset  of  these  points  in  constructing  an  actual  sample. 

Such  a  selection  night  be  made  in  many  possible  ways.  The  simplest  choice, 
and  the  only  one  we  have  investigated  to  date,  is  to  take  a  set  of  points 
uniformly  spaced  along  the  closed  search  curve.  We  require  that  one  of 
these  points  lies  at  s  =  0,  since  s  =  0  is  that  point  for  which  all  para¬ 
meters  k  take  on  their  nominal  values  kj[0'. 
i  t 

The  minimum  number  of  points  we  must  take  in  our  sample  can  be  related 

to  the  maximum  frequency  u)max  of  the  frequency  set.  This  relation  can  be 

derived  by  appeal  to  Nyquist’s  criterion  for  the  evaluation  of  Fourier 
(151 

coefficients'  .  We  want  to  evaluate  all  complex  Fourier  coefficients 
for  the  frequency  range 

°‘bl*2*w  •  (3.6) 

There  are  a  total  of  4w  +  1  such  coefficients.  By  Nyquist's  criterion, 

Ilia  X 

at  least  4ui  „  +  1  points  must  be  taken  on  the  search  curve  in  order  to  be 
max 

able  to  evaluate  this  number  of  coefficients. 
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It  also  is  convenient  to  let  the  number  of  points  be  even.  We  there¬ 


fore  let  N  =  2r  (>  4w  +  2)  be  the  number  of  uniformly  spaced  sample 

max 

points  along  the  search  curve,  and  define 

Sj  =  (j  =  1,2, ...2r)  .  (3.7) 

Furthermore,  if  we  choose  the  frequencies  w.  to  be  odd  integers  it  is 
easy  to  show  (Appendix  5)  that  the  output  function,  as  a  function  of  s, 
exhibits  the  symmetries 


fU/2  +  s)  =  f (tt/2  -  s) 
f (-ir/2  -  S)  =  f (-tt/2  +  S) 


(3.8) 


so  that  the  search  variable's  range  may  be  restricted  to  -ir/2  s  s  s  w/2 . 
Thus,  it  is  sufficient  to  evaluate  f(Sj)  only  for  those  s^  which  satisfy 


-ir/2  £  S.  £  tt/2.  (3.9) 

J 

The  number  of  points  which  satisfy  this  criterion  is  r  instead  of  2r.* 

If  we  choose  2r  to  be  of  the  form  4q  +  2  with  q  integer,  l.e.  divisible  by 
two  but  not  by  4,  then  eq.  (3.7)  can  be  modified  to 


or,  more  conveniently, 


r+1  r+3 

T~’  T"’ 


3r-3 


3r-l 


(3.10) 


s 


i 


it  /2z-r-l 
r  {~1 — 


■)  ; 


i 


1,  2, 


r 


(3.11) 


*  In  references  (1),  (2),  and  (3),  this  symmetry  was  not  invoked,  so  that 
the  number  of  sample  points  listed  in  these  papers  are  all  too  large  by 
a  factor  of  two. 
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Aside  from  a  lower  limiting  value*  r  a  2a,rn_v  +  1,  fixed  by  the  Nyquist 
cr  iterion,  r  can  be  assigned  any  value.  However,  large  values  of  r  are 
numerically  desirable  for  reasons  of  accuracy,  although  smaller  values 
are  desirable  for  reasons  of  computing  economy.  To  date,  in  practical 
applications  we  have  tended  to  allow  considerations  of  economy  to  prevail, 
and  we  usually  have  chosen  r  =  2w  +1. 

It  is  important  to  realize  that  a  more  accurate  evaluation  of  the 
Fourier  coefficients,  made  possible  by  choosing  r  to  be  large,  may  not 
result  in  improved  accuracy  for  the  sensitivity  coefficients.  For,  if 
the  search  curve  does  not  do  an  adequate  job  of  covering  the  parameter 
space,  then  increasing  the  accuracy  of  the  s-space  integration  will  not 
improve  the  coverage  of  the  parameter  space.  Thus,  in  addition  to  insuring 
that  r  is  sufficiently  large  to  obtain  an  accurate  s-space  quadrature,  we 
must  also  investigate  how  to  obtain  an  s-space  curve  which  does  a  good  job 
of  covering  the  entire  parameter  space. 

3.  Selection  of  Frequency  Set 

An  adequate  coverage  of  the  parameter  space  is  ensured  by  an  appropriate 
choice  of  the  integer  frequencies  and  the  number  of  points  2r.  Criteria 
for  such  a  choice  were  discussed  in  references  (1)  and  (2)  and  sets  of  u^,N 
(N=2r)  generated  according  to  these  criteria  are  presented  there.  Note 
that  once  a  set  nas  been  constructed  for  n  parameters,  this  same  set  can  be 
used  in  all  n  parameter  problems;  i.e.,  frequency  sets  can  be  tabulated 
once  and  for  al 1 . 

*  In  reference  (1),  it  ^s  indicated  that  by  appropriately  choosing  the 

frequency  set  it  s  possible  to  use  a  slightly  smaller  lower  limit 

value  for  r,  r*2<u  -  7.  The  practical  difference  between  these  two 

max 

limits  is  small,  and  we  prefer  to  use  the  larger  value  for  purposes  of 
discussion.  In  actual  calculations,  we  have  used  both  values. 


For  a  given  choice  of  frequency  set  w  and  number  of  points  2 r  we 
obtain  an  approximate  sensitivity  coefficient  which  we  designate  by  an 
asterisk.  Thus,  we  can  write 


A*  «  Aa  +  e^u^r)  *  (3.12) 

where  A  is  the  true  Fourier  coefficient,  A*  the  calculated  Fourier 
i  i 

coefficient  and  e  is  the  error,  which  depends  on  u,  2r  and  the  index  i  of 
the  coefficient  to  be  calculated. 

The  analysis  of  reference  (3)  and  parts  of  section  2  of  this  paper 
address  the  difficult  problem  of  how  to  determine  what  the  error  is  under 
given  circumstances  and  how  to  minimize  this  error  by  properly  choosing 
2r  and  w.  For  the  purpose  of  presentation  of  this  implementation,  we 
assume  that  the  error  has  been  made  sufficiently  small  to  carry  out  a 
meaningful  sensitivity  analysis. 

4.  Working  Equations 

The  symmetries  exhibited  in  eq.  (3.8)  simplify  the  computation  of  the 
numerical  Fourier  coefficients.  If  we  recall  that  the  total  number  of 
points  2r  is  chosen  as  4q+2  with  q  an  integer,  the  use  of  these  symmetries 
leads  to  (see  /Appendix  VI) 


A*«0 

i 

B*  =  0 
l 

AJ  *  (2q+l)_1 
B*  *  (2q+l)'1 


( i  odd) 
(i  even) 


U  odd) 


U  even) 


(3.13a) 

(3.13b) 

(3.13c) 


(3.13d) 
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where  we  have  set  f(s.)  =  fj.  These  formulae  generate  q+1  unique  cosine 

J  J 

coefficients  and  q  unique  sine  coefficients.  They  form  the  working  formulae 

of  our  method,  whereby  we  use  the  sample  values  f .  for  -q  s  j  s  q  to 

J 

generate  the  Fourier  coefficients. 

We  might  particularly  note  that  Eq.  (3.13c)  gives,  for  i  =  0 

q 

A*=(2q+1)'1  >  fi  (3.14) 

j=-q  J 

identifying  Ag  as  the  mean  value  of  f(s).  The  variance  of  f  is  generated 
by  the  formula 

q 

0*  =  (2q+l  )”^  l  (frC)Z  (3.15) 

j  =  -q  J 

It  is  worth  pointing  out  that  in  any  practical  calculations,  the 
amount  of  computing  needed  for  evaluating  Fourier  coefficients  is  a  very 
small  part  of  the  total  computing.  No  significant  benefit  is  thus  gained 
by  attempting  to  introduce  methods  such  as  Fast  Fourier  Transforms  (FFT) 
into  the  analysis.  There  would  in  fact  appear  to  be  disadvantages  in 
using  FFT  in  this  application,  since  FFT  works  best  when  the  numbers  of 
sample  points  is  highly  composite,  e.g.,  a  number  of  the  form  of  2n,  whose 
n  is  some  integer,  whereas  our  analysis  gains  in  simplicity  when  the  number 
of  points  is  not  highly  composite. 


5.  A  Useful  Modification 

Under  certain  circumstances,  partial  information  may  be  available 
concerning  the  dependence  of  an  output  function  upon  one  or  more  of  the 
parameters.  It  may  be  possible  to  use  this  information  to  improve  the 
numerical  accuracy  of  the  computation  of  the  partial  variances.  By  way 
of  example,  suppose  we  knew  that  a  particular  output  function  f  was  an 
even  function  of  a  particular  parameter,  say  ,  i.e.,  f(-k^)  *  f (k1 ) .  In 
such  a  case  it  is  easy  to  establish  that  the  Fourier  fundamental  of  fre¬ 
quency  u> i ,  and  all  odd  harmonics,  vanish  identically.  However,  even  harmonics 
would  not  in  general  vanish.  If  we  were  to  compute  on  the  basis  discussed 
previously,  half  of  the  terms  contributing  to  the  partial  sensitivity  would 
vanish,  including  the  fundamental  --  for  which  we  might  otherwise  expect 
higher  numerical  accuracy  than  for  the  even  harmonics. 

To  circumvent  this  loss  of  accuracy,  we  can  utilize  the  following 
trick.  We  define  a  new  parameter  kf  by 

k{  =  k*  (3.16) 


with  a  distribution  function  Pj’(kf)  related  to  P^(k^)  by 

dk,  P.Mc?) 

P{(k{)  =  P^k,)  -1=  -I— L 

1  i  i  1  ji.'  o  nr? 


(3.17) 


We  then  utilize  k|  and  Pj(k^),  instead  of  k^  and  P^(k^),  to  determine 
by  Eq.  (2.6).  Finally,  k^'  instead  of  k^  is  related  to  s  vi£  Eq.  (3.1). 

With  this  modification,  the  fundamental  Fourier  coefficient  ,  and 
in  general  all  of  the  harmonic  coefficients  do  not  vanish,  and  we  can 
expect  a  more  accurate  evaluation  of  the  partial  variance. 
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In  a  more  general  vein,  suppose  we  have  a  priori  knowledge,  or  suppose 
we  have  qualitative  reason  to  believe,  that  the  output  function  f  depends 
upon  k-j  via  an  explicit  function  of  k^ .  This  is  to  say,  suppose  f(k^)  is 
actually  of  the  form  F(h(k^)),  with  h(k-j)  explicit.  Then  it  is  best  to 
proceed  by  defining 


Pf(kf)  =  P-j  (k-| ) 


kj-  =  h (k-j ) 


dk1  P1(k1)  P1 

h'hk,)] 

h  '(k1 )  h. 

1 

*■ 
— *  \ 

1  .  J 

(3.18) 


Next  utilize  P,'(kT)  to  determine  G-|  by  Eq.  (2.6),  and  relate  k-j”  to  s  by 
Eq.  (3.1).  After  this,  proceed  along  the  lines  previously  discussed. 
The  use  of  these  modifications  can  be  expected  to  improve  the  numerical 
accuracy  of  the  sensitivity  analysis. 
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IV.  Appl i cations 

To  illustrate  the  application  of  our  sensitivity  analysis  method,  we 
present  in  this  section  four  different  examples.  The  first  three  of  these 
examples  have  not  been  published  previously,  although  two  of  them  are 
available  through  government  document  sources. 

The  applications  to  be- discussed  are  the  following: 

1.  An  Economic  Model 

2.  A  Chemical  Laser  Model 

3.  A  Socioeconomic  Model  (World  II) 

4.  A  Chemical  Reaction  Model 

In  the  first  of  these  applications  we  are  dealing  with  a  model  which  uses 
linear  programming  equations.  The  second  and  fourth  examples  deal  with 
models  described  by  ordinary  differential  equations.  The  third  application 
involves  a  set  of  differential  difference  equations.  Our  ability  to 
achieve  a  successful  sensitivity  analysis  for  such  different  types  of  model 
equations  demonstrates  the  wide  applicability  of  our  technique  of  sensitivity 
analysis. 

One  important  feature  should  be  noted  which  pervades  the  analysis  of 
all  these  examples.  A  study  of  the  sensitivity  coefficients  reveals  un¬ 
expected  but  significant  relations  between  parameters  and  output  functions 
which  could  not  have  been  predicted  from  a  more  conventional  analysis  or  a 
purely  intuitive  approach.  Sensitivity  analysis  thus  provides  information 
which  can  lead  to  important  insights  into  the  structure  of  the  models 
used  to  represent  complex  systems. 

1.  AN  ECONOMIC  MODEL^ 

Frequent  use  Is  made  of  linear  programming  (LP)  and  related  methodologies 
in  modeling  the  operation  of  complex,  interacting  systems.  Generally 
speaking,  the  form  in  which  such  problems  are  posed  is  one  of  seeking  a 
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constrained  optimum  to  some  objective  function.  Often  the 
'■‘t  objective  function  is  the  profit  to  be  derived  from  a  business, 
and  the  constraints  relate  to  limitations  on  available  re¬ 
sources,  manpower,  capital,  equipment,  and/or  demand.  There 

is  an  enormous  literature  covering  this  subject,  of  which  we 

[7] 

only  cite  a  few  examples.  J  An  interesting  particular 
example  of  application  to  a  model  of  operation  of  a  petroleum 
refinery  is  described  in  a  text  by  Wilde  and  Beightler. 1  x  The 
structure  of  this  model  is  such  that  it  can  be  extended  to 

characterize  an  aggregated  representation  of  the  entire  domestic 
petroleum  industry.  In  this  model,  a  number  of  crude  petroleum 
producers  supply  crude  oils,  each  of  its  own  characteristic 
chemical  composition,  to  a  number  of  refineries,  each  of 
specific  design.  The  different  crudes  each  have  a  range  of 
possible  products  into  which  they  may  be  converted,  and  each 
refinery  has  a  range  of  operability  with  respect  to  these 
crudes.  For  each  crude,  the  supply  is  limited  in  an  absolute 
sense  by  the  pumping  capacity  of  the  supplier.  For  each  re¬ 
finery,  the  process  capacity  is  limited  by  existing  equipment 
limitations.  In  addition  to  these  physical  constraints,  there 
are  economic  constraints.  The  crude  suppliers  may  limit  pro¬ 
duction  to  below  their  physical  capacity  in  order  to  obtain 
political  or  economic  advantages.  The  processors  may  limit 
production  for  similar  reasons.  Furthermore,  there  may  be 
production  limits  due  to  consumers,  who  at  any  given  time 
have  limitations  upon  their  desire  to  use  or  capacity  to  store 
products . 
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All  of  these  factors  and  others  serve  to  de¬ 
termine  the  rates  of  consumption  of  various  crude  petroleums, 
the  rates  of  production  of  various  products,  and  the  rates  of 
consumption  of  these  products.  The  money  flow  from  consumer 
to  producer  to  supplier  is  likewise  determined  by  these  factors, 
as  are  the  profits  of  the  producers  and  suppliers.  Associated 
with  many  of  these  factors  are  numerical  values  of  parameters, 
almost  all  of  which  are  imprecisely  known.  Furthermore,  these 
parameters  are  not  static  but  tend  to  fluctuate  in  time,  both 
in  response  to  factors  internal  to  the  petroleum  industry,  and 
to  socioeconomic,  political,  and  technical  factors  external 
to  the  petroleum  industry.  These  fluctuations  make  profit 
optimization  difficult,  and  make  performance  prediction  for 
the  industry  equally  difficult. 

The  objectives  of  our  study^  were  to  develop  a  simple  aggregated 
model  of  industry  economics,  apply  sensitivity  analysis  to  the  model,  and  to 
see  whether  the  resulting  methodology  might  be  useful  as  a  means  of  determining 
which  parameters  most  influence  costs,  profits,  and  the  development  of  shortages. 

The  mathematical  structure  of  the  model  is  as  follows: 

The  symbol  xiJc  denotes  the  daily  national  utilization  of  crude 
petroleum  available  from  source  i  and  which  is  converted  to 
products  by  process  k.  The  physical  constraints  of  supply 
can  be  written  as 
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< 


£*ik<bi  <i-l.  2 . N)  (4.1) 

k 

where  is  the  availability  of  type  i  crude.  The  constraints 
of  consumer  demand  are  of  the  type  *• 

N 

yi  S  ajik  Xik  —  faN+j 
i=l  k 


(j  -  l,  2,  ...,  Mi  ;  (4.2) 


The  meaning  of  this  relation  is  that  production  will  not  exceed 
demand.  Clearly  this  relation  is  not  instantaneously  true,  but 
rather  represents  a  time  average.  That  is,  at  any  given  instant, 
production  may  well  exceed  demand,  but  over  the  long  term,  either 
demand  will  increase  (perhaps  due  to  price  adjustments)  or  pro¬ 
duction  will  be  reduced  so  as  to  restore  a  condition  in  which 
the  inequality  (4^.2)  is  satisfied. 
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where  ajj_jc  the  amount  of  product  j  produced  from  crude  i 
by  process  k,  and  bN+j  is  the  consumer  demand  for  product  j. 

The  production  industry  tries  to  operate  so  as  to  opti¬ 
mize  the  daily*  profit  P,  which  is  given  by 


M 


P  ■  12  Vj  ■  12  xipi 

j-i 


N 

E 

i=l 


(4.3) 


where 


i-E 


xik 


(4.4) 


and  where  s^  is  the  selling  price  of  product  j,  y^  is  the  daily 
production  of  product  j,  x^  is  the  daily  production  of  crude  i, 
and  p^^  is  the  cost  of  crude  i.  The  latter  cost  includes  the 
cost  of  transportation  to  the  refinery,  and  an  allocated  cost 
of  shipping  to  the  consumer  for  each  product  manufactured  from 
crude  i  at  the  refinery.  Production  and  utilization  are  re¬ 
lated  by 


N 


EE 

i=k  k 


ajik  xik 


(*.5) 


so  that  by  substitution 

*  ■  1212  f'pi + 12  *j ik 

i  k  1  j 


ik 


(4.6) 


We  use  in  general  one  day  as  the  unit  of  time,  one  barrel  as 
the  unit  of  production,  both  of  crudes  and  of  products,  and 
one  dollar  as  the  unit  of  profit  or  loss. 
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For  a  given  set  of  prices  consumer  demand  will  vary. 
Although  the  finer  details  of  this  variation  may  be  obscure, 
it  is  plausible  to  model  this  variation  by  the  relations 


where  g^  is  the  limiting  consumer  demand  as  the  price  falls  to 

zero,  and  s^u  is  the  limiting  upper  price  at  which  demand 

* 

falls  to  zero.  The  exponent  fixes  the  marginal  demand, 
i.e.,  the  rate  of  change  of  demand  with  respect  to  price. 

Small  values  of  are  referred  to  as  indicating  large 
elasticity  of  demand,  and  large  values  of  q^  are  referred  to 
as  indicating  small  elasticity  of  demand.  For  essential  goods, 
such  as  many  petroleum  products,  q^  tends  to  be  large  (i.e., 

>>  1).  The  specific  values  of  g ^ ,  s^u  and  q^  for  a  given  pro¬ 
duct  must  be  regarded  as  experimental  parameters  which  can 
only  be  determined  by  observations  of  the  response  of  the 
market  place  to  price  changes.  In  addition,  they  are  not  truly 
constants  for  any  product,  but  display  variations  in  time  due 
to  changes  in  technology,  social  structure,  political  institu¬ 
tions,  public  policies  and  other  factors. 

If  we  substitute  Eq.  ( 4jL  7. )  into  Eq.  (4J2),  we  are  then 
confronted  with  the  problem  of  maximizing  P  as  a  function  of 
the  prices  s^  and  the  production  levels  x^,  subject  to  the 
* 

Equation  (4.7)  is  an  approximate  form  which  should  not  be  in¬ 
terpreted  as  literally  valid  over  too  broad  a  range  of  s^. 

If  product  i  is  an  essential  material,  nations  will  go  to 
war  at  prices  well  below  Sj.  *  s^u,  if  war  promises  to  induce 
a  cheaper  supply.  Also,  even  if  s^  were  to  approach  zero, 
demand  does  not  become  infinite  because  of  limits  in  con¬ 
sumer's  storage  capacity. 


*43- 


nonlinear  constraints,  Eq.  (4.2),  and  the  linear  constraints 
* 

Eq.  (  4.1).  This  problem  is  nonlinear  and  can  be  solved.  How¬ 
ever,  we  choose  instead  to  introduce  a  simplification  which  re¬ 
duces  the  problem  to  a  linear  one.  To  do  this,  we  require 
that  all  products  be  profitable,  a  condition  we  assure  by  re¬ 
quiring  all  product  prices  to  equal  or  exceed  the  price  of  the 
most  expensive  crude,  with  a  markup  to  account  for  processing 
losses.  To  state  this  mathematically,  we  define 


j 


that  is,  a^  is  the  fraction  of  each  barrel  of  crude  i  lost 
during  production  by  process  k.  Then  define  a  lower  price 
level  Sj  by 

s «  *  ““  — ^ -  (j  =  1,  2,  ...,  N)  (4.9) 

'  U-aik> 

Sale  of  any  product  at  a  price  (per  barrel)  in  excess  of  sA 
will  be  profitable. 

Define  then  a  parameter  A^  for  each  product  such  that 
its  price  is  given  by 

si  ”  si  +  Xi<siu~s!}*  2>  •••#**)  (.4.10) 


Then  if  A^.  *  0,  we  have  s^^  *  s4?  but  if  A^^  *  1,  we  have 
sA  ■  s^u,  the  limiting  upper  price  at  which  sales  of  product 
i  would  fall  to  zero. 


This  statement  views  the  problem  from  the  point  of  view  of  the 
refiners.  In  reality,  the  problem  is  more  complex,  because 
the  producers  interest  is  in  maximizing  his  profit  (which  is 
not  expressed  by  Eq.  (4.6))  subject  to  the  constraints.  The 
producer's  tool  for  attempting  this  is  by  variation  of  the 
crude  prices  pj,  analysis  of  which  presents  a  more  complex 
problem  than  tnat  we  will  consider  here. 
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* 

Clearly  then  each  X^  must  fall  into  the  range 
0  <  X  <  1.  We  can  refer  to  the  X^'s  as  the  "markup  para¬ 
meters"  for  the  products,  and  utilize  the  X^'s  as  equivalent 
to  a  representation  of  the  prices  s^. 

We  can  develop  an  analysis  in  which  the  X^'s  are  varied 
independently  of  the  crude  prices  p^ .  If  we  then  consider  a 
fixed  set  of  X^'s,  the  prices  and  hence  the  consumer  demand 
constraints,  Eq.  (4.2),  will  be  determined  via  Eqs.  (4~I7)  - 
(47.10).  .Then  for  a  fixed  set  of  supply  constraints,  Eq.  (4.J1), 
the  profit  optimization  problem  reduces  to  a  linear  program¬ 
ming  problem,  namely  that  of  maximizing  P  subject  to  the  linear 
constraints  of  Eqs.  (4.1)  and  (4.2).  The  form  of  this  problem 
is  essentially  the  same  as  that  discussed  by  Wilde  and 
Beightler. ^ 

We  now  have  defined  a  linear  programming  problem  in 
which  producer's  seek  to  maximize  the  total  profit  P,  Eq.  (4~.6) 
for  given  values  of  crude  prices  p ^ ,  crude  availabilities  b^ , 
product  saturation  demands  g^,  product  cutoff  prices 
product  markup  rates  X^  and  demand  exponents  q^.  None  of  these 
parameters  are  constants,  but  rather  they  all  vary  from  time 
to  time  in  response  to  social,  economic  and  psychological 
forces.  It  is  important  to  try  to  understand  how  the  variabil¬ 
ity  of  supply  and  price  will  influence  costs,  consumptions, 

. 

"It  Is  Implicit  that  s^  s  n\n  s^u,  which  effectively  states  that  we  consider 
only  the  economic  circumstances  in  which  producer  prices  are  not  set  so  high 
as  to  prevent  sales  of  an^  product  because  of  the  impact  of  Eq.  (4.7). 
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and  shortages  of  refined  products.  Markup  rates  are  subject 
to  fluctuation.  Saturation  demands,  cutoff  prices,  and  demand 
exponents  are  only  measurable  approximately,  and  furthermore 
fluctuate  in  time.  Thus  we  have  an  LP  problem  with  a  large 
number  of  parameters  of  uncertain  value.  The  problem  solu¬ 
tion  depends  upon  these  parameter  values.  What  we  have  done 
is  to  treat  these  parameter  uncertainties  statistically,  there¬ 
by  generating  a  set  of  linear  programming  problems.  We  have 
performed  a  sensitivity  analysis  and  determined  which  para¬ 
meter  uncertainties  are  most  influential  in  causing  varia¬ 
bility  in  the  following  dependent  variables:  product  con¬ 
sumptions,  crude  consumptions,  expenditures  for  crudes,  ex¬ 
penditures  for  products,  profits  and  unfilled  demands,  i.e., 
shortages . 

In  the  study  that  was  performed,  only  a  subset  of  the 
parameters  was  taken  to  be  uncertain:  crude  prices,  product 
markup  rates  and  the  demand  exponent.  It  was  furthermore 
assumed  that  the  demand  exponent  was  the  same  for  all  pro¬ 
ducts,  i.e., 

qi  *  q  (i  *  1,  ...,  M)  14.11) 

Other  parameters  were  given  constant  values  (i.e.,  no  van- 
certainty  was  assumed)  based  upon  a  brief  analysis  of  data 
on  petroleum  economics.  The  range  of  prices  was  such  that 
price  and  not  physical  supply  constrained  production. 

To  illustrate  application  of  the  methodology,  an  ag¬ 
gregated  model  was  chosen  in  which  there  are  four  types  of 
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crudes  available.  Of  these,  three  could  be  processed  only 
in  one  way,  and  the  fourth  had  two  alternative  processes 
applicable  to  it.  The  crude  production  variables  then  are 

labeled  x^,  x2i'  x3i/  x4i/  and  x42*  Pro<^uc'ts  were  put  into 
four  classifications,  termed  "gasoline"  (subscript  *=1), 
heating  oil  (subscript  =  2),  lubricating  oil  (subscript  *  3), 
and  aviation  fuel  (subscript  =  4) .  Available  crude  supplies, 
product  saturation  demands  and  product  cutoff  prices  are 


given  in  Table  4~.  1 . 
Table  4.2. 


The  coefficient  matrix  a 


jik 


is  given  in 


The  ranges  of  the  parameter  values  are  shown  in 
Table  4.3.  Because  we  did  not  have  statistical  data  on  the 


probability  distribution  of  these  parameters,  we  in  all  cases 
assumed  a  uniform  distribution  between  the  indicated  lower 
and  upper  limits.  All  of  the  parameter  ranges  were  based 
upon  rough  estimates  derived  from  real  data.  The  crude  prices 
quoted  in  Table  4.3  includes  costs  of  transportation  to  the 
refinery,  costs  of  product  transportation  to  the  retail 
market,  and  costs  of  processing.  The  sum  of  these  costs  is 
about  $5. 00/barrel,  so  that  the  range  of  wellhead  prices  which 
corresponds  to  Table  4.3  is  about  $7.00  to  $11. 00/barrel 
($5.00  to  $10. 00/barrel  for  Type  3  crude). 

With  nine  uncertain  parameters,  it  is  necessary  to  solve 
323*  linear  program  problems,  using  323  parameter  sets,  as 

*According  to  reference  2,  630  problems  must  be  solved,  but  the 

two-fold  symmetry  discussed  in  Section  Illreduces  this  to  315 

problems.  The  slightly  higher  value  323  is  obtained  from  the 

relation  N  »  21T  +1  *  2  (161)  +  1  *  323. 

max 
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TABLE  4.1 

CONSTRAINTS  AND  PARAMETERS  OF  A  TEST  PROBLEM  IN 
PETROLEUM  ECONOMICS 

Available  Crude  Supplies 

='7.0  *  10®  barrel/day 

b2  =  7.0  x  io®  barrel/day 

b^  =  7.0  x  io®  barrel/day 

b4  =  4.2  x  io®  barrel/day 

Product  Saturation  Demands 

g^  =  9.8  x  10®  barrel/day 

g2  =  4.9  x  io®  barrel/day 

=  6.0  x  io5  barrel/day 

g4  =  2.1  x  io®  barrel/day 

s^u  =  $42. 00/barrel 
s2u  *  $46 .00/barrel 
s3u  =  $180 . 00/barrel 
s4u  =  $60 . 00/barrel 


TABLE  4.3 

RANGE  OF  VALUES  OF  UNCERTAIN  PARAMETERS 


Crude  Prices 

Pi 

P3 

P4 


$12 . 00/barrel 
$12. 00/barrel 
$10. 00/barrel 
$12. 00/barrel 


to  $17 .00/barrel 
to  $17 . 00/barrel 
to  $15. 00/barrel 
to  $17 .00/barrel 


Product  Markup  Parameters 


0.03  to  0.07 
\2  0.03  to  0.07 
X3  0.03  to  0.07 
X^  0.02  to  0.04 


Product  Markup  Exponent 


q  5.9  to  6.1 


I 


, 

defined  by  Eq.  (3.1).  The  nine  frequencies  are  given  in 

Table  I  of  reference  (2). A  computer  program  was  written  which 

selects  parameter  sets  and  which  then  calls  on  a  standard  linear 

programming  subroutine  to  solve  a  problem.  The  results  are 

then  collected  and  transformed  by  the  main  program  so  as  to 

2  - 

provide  mean  values,  <f>.  variances  and  partial  variances  j 
S*  (c.f.)Eqs.  .(3.2)and{3.15)flh®  output  functions  studied  were  the 

l 

consumption  rates  for  the  individual  crude  petroleums,  the 
daily  expenditures  for  all  crudes,  and  for  all  products,  the 
daily  profits  of  the  refineries,  the  daily  deliveries  of  the 
products  and  the  unfilled  demands.  Tables  4.4,  4.5,  and  4.6 
summarize  the  results. 


I 


In  Table  4.4  we  consider  the  variables  x^  (c.f.,  Eq. 
(4\4)),  the  consumptions  of  the  four  crude  petroleums.  For 
each  crude,  the  table  lists  the  available  3 

supply.  Also  listed  are  the  minimum  and  maximum  consumption 
as  obtained  from  the  sample  of  323  calculations.  In  all 
cases  these  vary  from  the  maximum  available  down  to  zero. 

The  nominal  consumptions  correspond  to  the  single  choice 
s  »  0  for  the  search  variable.  For  a  uniform  probability  dis¬ 
tribution  this  means  that  all  variables  have  been  taken  as  the  average 
value  in  their  ranges  of  uncertainty  (c.f..  Table  4.3).  The  coefficient  of 
variation  v^  provides  an  alternative  representation  of  the  variance.  By 
definition  the  coefficient  of  Variation,  v^,  of  a  random  variable  f  is  given 
by 

vf  *  of/<f>  (4.12) 
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TABLE  4.4 


SENSITIVITY  ANALYSIS  OF  CONSUMPTIONS  OF  CRUDE  OILS, 
DAILY  EXPENDITURES,  AND  PROFIT  RATES 


Crude  Consumptions  (bbl/day) 


Crude  Type  1 

Crude  Type  2 

Crude  Type  3 

Crude  Type 

Available  Supply 

7.0  x  106 

7.0  x  106 

7.0  x  106 

4.2  x  106 

Minimum  Daily  Consumption 

0 

0 

0 

0 

Nominal  Daily  Consumption 

7.0  x  106 

1.3  x  106 

7.0  x  106 

3.0  x  106 

Average  Daily  Consumption 

4.7  x  106 

3.3  x  106 

6.1  x  106 

3.6  x  106 

Maximum  Daily  Consumption 

7.0  x  106 

7.0  x  106 

7.0  x  106 

4.2  x  106 

Coefficient  of  Variation 

0.686 

0.921 

0.249 

0.202 

Partial  Variances  $* 

“o 

A, 

Price  of  Crude  No.  1  (p^) 

0.422 

0.124 

0.025 

0.012 

Price  of  Crude  NO.  2  (p2) 

0.197 

0.567 

0.275 

0.247 

Price  of  Crude  No.  3  (p3) 

0.022 

0.039 

0.252 

0.001 

Price  of  Crude  No.  4  (p^) 

0.004 

0.011 

0.001 

0.398 

Markup  of  Product  1  (A^ 

0.004 

0.000 

0.003 

0.006 

Markup  of  Product  2  (Aj) 

0.001 

0.000 

0.000 

0.001 

Markup  of  Product  3  (A^) 

0.004 

0.001 

0.000 

0.002 

Markup  of  Product  4  (A4) 

0.004 

0.001 

0.001 

0.002 

Product  Markup  Exponent  (q) 

0.003 

0.001 

0.006 

0.003 
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The  data  in  Table  4.4  indicates  considerable  variability  in 
consumption  of  the  individual  crudes  as  a  consequence  of 
changes  in  the  parameters.  The  partial  variances  indicate 
the  differing  degree  to  which  parameter  uncertainties  con¬ 
tribute  to  consumption  variability.  We  briefly  summarize  the 
much  longer  discussion  of  results  contained  in  the  report  by 
Levine. ^  First,  it  is  not  surprising  that  consumption  of  Crude 
Type  1  is  more  sensitive  to  variation  in  the  price  of  Crude 
Type  1  than  it  is  in  any  other  variable.  The  only  other 
variable  of  importance  in  this  regard  is  the  price  of  Crude 
Type  2.  Similarly,  consumption  of  Crude  Type  2  is  most  sensi¬ 
tive  to  the  cost  of  Crude  Type  2.  However,  when  we  come  to 
Types  3  and  4  Crude  the  situation  is  different.  Consumption 
of  Type  3  is  more  dependent  upon  price  of  Type  2  than  on  its 
own  price,  arid  consumption  of  Type  4  is  nearly  as  sensitive 
to  the  price  of  Type  2  as  it  is  to  its  own  price.  These 
latter  relations  suggest  a  tendency  for  Type  2  to  serve  as  a 
substitute  for  Types  3  and  4  under  certain  conditions  of 
relative  price. 

In  Table  4.5  we  consider  costs  of  crudes  and  of  pro¬ 
ducts,  and  refinery  profits.  The  table  clearly  indicates 
that  the  net  expenditures  for  all  crude  petroleum  is  more 
sensitive  to  the  price  of  Type  3  Crude  than  to  any  other 
parameters,  somewhat  less  sensitive  to  the  price  of  Type  2 
Crude,  and  comparatively  insensitive  to  the  other  parameters. 

This  conclusion  with  respect  to  Type  3  Crude  is  in  accord 


TABLE  4.5 


SENSITIVITY  ANALYSIS  OF  EXPENDITURE  RATES  FOR 
CRUDE  PETROLEUMS,  EXPENDITURE  RATES  FOR 
REFINED  PRODUCTS,  AND  REFINERY  PROFITS 
(Units:  Dollars/Day) 

Expenditures  for  Expenditures  for 

Crude  Petroleum  Refined  Products  Profit 


Minimum 

$1,583 

X 

108/aay 

$2,301  x  108/day 

$2,864 

x  107/day 

Nominal 

2.511 

X 

108 

2.915  x  108 

4.046 

x  107 

Average 

2.358 

X 

108 

3.027  x  108 

6.695 

x  107 

Maximum 

2.827 

x  108 

8 

3.488  x  10 

1.119 

x  108 

Coefficient  of  Variation  v 

0.098 

0.076 

0.243 

Partial  Variances  S* 

JL 

Price  of  Crude  No.  1  . 

0.070 

0.043 

0.071 

Price  of  Crude  No.  2 

0.169 

0.328 

0.112 

Price  of  Crude  NO.  3 

0.375 

0.059 

0.287 

Price  of  Crude  NO.  4 

0.065 

0.133 

0.059 

Markup  of  Product  1  (X^) 

0.006 

0.040 

0.030 

Markup  of  Product  2  (X2) 

0.001 

0.009 

0.010 

Markup  of  Product  3  (X3) 

0.007 

0.012 

0.003 

Markup  of  Product  4  (X^) 

0.002 

0.002 

0.000 

Product  Markup  Exponent  (q)  0.006 


0.006 


0.001 


with  the  data  in  Table  4.4,  which  shows  the  daily  average 


consumption  of  Type  3  Crude  to  be  higher  than  that  of  any 

/ 

other  crude,  either  on  an  absolute  or  a  relative  basis.  Less 
obvious  is  the  ranking  of  the  price  of  Type  2  Crude  as  the 
second  most  important  parameter,  particularly  when  one  notes 
that  its  daily  average  consumption  is  the  least  of  all  of  the 
crudes.  We  also  note  from  Table  4  .4  that  the  coefficient  of 
variation  is  higher  for  the  consumption  of  Type  2  than  it  is 
for  any  other  crude,  so  that  its  average  consumption  repre¬ 
sents  the  mean  value  of  a  widely  dispersed  variable. 

The  observation  of  the  importance  of  the  price  of 
Type  2  Crude  as  a  parameter  is  a  good  illustration  of  the 
ability  of  sensitivity  analysis,  in  the  form  we  have  developed, 
to  locate  obscure  but  sicnificant  relationships.  By  con¬ 
trast,  if  we  had  only  considered  nominal  parameter  values,  or 
small  variations  about  nominal  parameter  values,  we  would 
have  found  the  nominal  daily  consumption  of  Type  2  Crude  to 
be  the  least  of  all  the  crudes  (see  Table  4  .4),  and  we  likely 
would  have  concluded  that  its  price  did  not  importantly  af¬ 
fect  net  expenditures.  Sensitivity  analysis  thus  gives  us 
clues  which,  if  pursued,  can  lead  to  deeper  insights  into 
the  structure  of  a  complex  system  than  we  could  obtain  by 
more  conventional  analysis. 

In  the  same  vein.  Table  4.5  shows  the  price  of  Type  2  Crude  to  be 
the  most  important  parameter  in  the  determination  of  expenditures  for 
refinery  products,  despite  the  fact  that  the  average  daily  consumption 
of  this  Crude  is  less  than  that  of  the  ether  Crudes. 


It  is  difficult  to  provide  a  simple  explanation  of  the 

high  importance  of  Type  2  Crude,  and  the  explanation  must 


lie  in  a  complex  interplay  of  the  relative  prices  of  the 
crudes,  the  different  product  distributions  available  from 
refining  different  crudes,  and  the  differing  product  demands. 
Nonetheless,  the  fact  of  this  high  importance  emerges  directly 
from  the  analysis. 

It  is  perhaps  worth  remarking  on  the  almost  complete 
insensitivity  of  the  dependent  variables  to  the  product  mark¬ 
ups  and  the  markup  exponent.  This  conclusion  is  in  accord 
with  the  real  world  observation  that  consumption  of  petroleum 
products  continues  unabated  even  in  the  face  of  large  price 
increases;  i.e.,  the  price-demand  relation  is  "highly  in¬ 
elastic",  in  the  language  of  economics.  Petroleum  products 
are  simply  too  important  to  the  consumer  for  their  use  to  be 
foregone  in  virtually  all  circumstances. 

Table  4.6  combines  our  sensitivity  analyses  for  product 
delivery  rates  and  for  shortages.  We  can  do  this  in  one  table 

because  shortage  is  defined  as  the  algebraic  difference  be- 

* 

tween  demand  and  delivery.  Because  this  relation  is  linear. 

In  constructing  Table  4.6,  we  can  define  "shortage"  in  two 
alternate  ways:  the  difference  betweel  actual  demand  at  the 
current  price  level  and  delivery;  or  by  the  difference  between 
saturation  demand  and  delivery.  In  practical  application  to 
products  with  an  inelastic  price  demand  curve,  these  two 
definitions  are  virtually  coincident.  This  is  the  situation 
that  applies  in  this  study.  For  the  sake  of  definiteness, 
we  define  shortages  relative  to  saturation  demand;  i.e., 
we  use  the  second  definition. 


1^"  1/  V'  -  - 


TABLE  4.5 


SENSITIVITY  ANALYSIS  OF  PRODUCT  DELIVERY  RATES 
AND  UNFILLED  DEMAND  (SHORTAGES) 

Product  Deliveries  and  Shortages  (Barrels) 


1 

2 

3 

4 

Gasoline 

Heating  Oil 

Lubricating  Oil 

Aviation  Fuel 

Product  Saturation  Demand  (Daily) 

9.8  x  106 

4.9  x  106 

6.0  x  105 

2.1  x  106 

.■Minimum  Daily  Delivery 

6.3  x  106 

4.1  x  106 

0 

1.4  x  106 

Minimum  Daily  Shortage 

2.0  x  104 

6.7  x  103 

0 

1.9  x  105 

Nominal  Daily  Delivery 

9.1  x  106 

4.9  x  106 

6.0  x  105 

1.8  x  106 

Ncminal  Daily  Shortage 

6.6  x  105 

1.5  x  104 

0 

2.7  x  105 

Average  Daily  Delivery 

8.7  x  106 

4.9  x  106 

5.9  x  105 

1.7  x  106 

Average  Daily  Shortage 

1.1  x  106 

4.2  x  104 

1.1  x  104 

4.0  x  105 

Maximum  Daily  Delivery 

9.8  x  106 

4.9  x  106 

6.0  x  10S 

1.9  x  106 

Maximum  Daily  Shortage 

3.6  x  106 

8.4  x  105 

6.0  x  105 

7.1  x  105 

Coefficients  of  Variation  v+ 

Daily  Delivery 

0.105 

0.022 

0.103 

0.071 

Daily  Shortage 

0.867 

2.539 

5.301 

0.303 

Partial  variances  S* 

JL 

Price  of  Crude  No.  1 

0.399 

0.031 

0.017 

0.373 

Price  of  Crude  No.  2 

0.076 

0.023 

0.072 

0.002 

Price  of  Crude  No.  3 

0.098 

0.064 

0.009 

0.048 

Price  of  Crude  No.  4 

0.000 

0.007 

0.079 

0.128 

Markup  of  Product  1  (^) 

0.008 

0.001 

0.002 

0.004 

Markup  of  Product  2  (*2) 

0.001 

0.003 

0.001 

0.003 

Markup  of  Product  3  (^) 

0.007 

0.008 

0.018 

0.009 

Markup  of  Product  4  (X^) 

0.004 

0.002 

0.000 

0.002 

Product  Markup  Exponent 

0.007 

0.015 

0.007 

0.008 
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the  partial  variances  of  the  delivery  for  a  particular  product 
will  be  the  same  as  the  partial  variances  for  the  shortages 
of  the  same  product. 

With  respect  to  aviation  fuel,  the  results  shown  in 
Table  4,6.  indicate  that  the  range  of  daily  deliveries  is  small, 
the  average  and  nominal  daily  deliveries  are  nearly  equal, 
and  the  coefficient  of  variation  is  small.  The  price  of 
Type  1  Crude  is  the  most  important  parameter  and  the  price 
of  Type-  4  Crude  is  the  next  most  important  parameter.  The 
rationale  for  the  latter  observation  is  difficult  to  establish: 
since  all  of  the  Crudes  can  be  used  to  produce  aviation  fuel 
in  the  same  relative  proportions,  we  might  expect  that  the  para¬ 
metric  importance  would  be  in  the  same  order  as  the  average  daily 
consumption,  but  comparison  of  Tables  4.4  ano^4.6  shows  this  not  to  be  the  cast- 
Thus  we  again  see  an  example  where  the  numerical  results  of  a  sensi¬ 
tivity  analysis  do  not  conform  to  those  derived  by  simple  reasoning 

With  respect  to  gasoline,  a  more  complex  situation 
arises:  there  is  a  moderately  large  spread  of  delivery  rates, 
and  a  corresponding  spread  of  shortages  which  yields  at  the 
upper  extreme  shortages  which  are  a  large  fraction  37  per¬ 
cent)  of  demand.  The  only  parameter  of  major  significance 
influencing  this  variability  is  the  price  of  Type  1  Crude, 
despite  the  fact  that  it  ranks  second  in  average  daily  con¬ 
sumption  . 

With  respect  to  heating  oil,  still  another  situation 
arises.  The  range  of  daily  delivery  is  small  (as  reflected 


1 


in  a  small  coefficient  of  variation) ,  but  the  coefficient  of 
variation  of  the  daily  shortage  is  large.  The  latter  fact 
stems  from  the  evaluation  of  shortage  as  the  small  difference 


between  two  large  numbers  ,  demand  and  delivery  .  On  an 
absolute  basis,  the  maximum  shortage  is  17  percent  of  satura¬ 
tion  demand,  but  the  average  shortage  is  less  than  1  percent 
of  saturation  demand.  For  this  case,  the  sensitivity  analysis 
does  not  clearly  identify  a  parameter  of  maximum  importance, 
and  the  sum  of  the  partial  variances  in  Table  4.6  is  only 
0.151,  much  less  than  the  maximum  possible  value  of  unity. 


This  means  that  the  higher  order  partial  variances  oj^ , 
etc.  (c.f.,  Eqs.  (2.28b)  and  (2.28c))  are  important  in  this 


case.  Although  we  could  calculate  these  higher  order  partial 


variances,  we  have  chosen  not  to  do  so,  since  the  relatively 


low  level  of  shortages  renders  immaterial  such  a  complicated 


analysis. 

The  situation  for  lubricating  oil  is  similar  to  that 
for  heating  oil,  except  that  the  range  of  daily  shortages  is 
even  broader,  being  at  the  one  extreme  zero  and  at  the  other 
equal  to  the  product  saturation  demand.  The  sum  of  the  par¬ 
tial  variance j  in  Table  4.6  is  0.205,  again  well  below  unity 
in  value.  Here,  evaluation  of  higher  order  partial  variances 

would  be  useful,  so  as  to  help  identify  those  parameter  combi- 

* 

nations  which  induce  large  shortages. 


*This  identification  can  also  be  made  by  direct  examination  of 
the  323  sample  calculations,  which  would  be  a  tedious  process. 


SJ 


2. 


A  CHEMICAL  LASER  MODEL 


— 

The  hydrogen  fluoride  chemical  laser  is  a  chemical 
system  consisting  initially  of  a  mixture  of  molecular  hydro¬ 
gen,  H2»  molecular  fluorine,  F^,  and  an  inert  diluent,  e.g., 
argon.  At  time  t  =  0  flashlamp  irradiation  dissociates 

some  of  the  F2  into  F  atoms.  Chain ^reactions  then  _ 1... 

gin  which  leads  to  the  conversion  of  the  hydrogen  and  fluorine 
into  hydrogen  fluoride,  HF.  Because  of  the  high  chemical 
energy  release  the  HF  molecules  produced  are  in  comparatively 
high  vibrational  quantum  states.  If  the  system  is  placed  be¬ 
tween  two  mirrors,  the  vibrational  decay  will  lead  to  lasing 
action  on  the  infrared  vibration-rotation  transitions  of  the 
HF  molecule.  Theoretical  descriptions  of  the  system  have 
been  provided  by  various  authors.  .12]  We  have  performed  a 

sensitivity  analysis  on  _  the  model  of  Kerber,  Emanuel, 
and  Whittier. ^ 

The  model  studied  described  the  chemical  evolution  of 
the  system  in  terms  of  a  set  of  68  reversible  chemical  re¬ 
actions  between  the  following  species:  H  atoms,  Ar  atoms, 

F  atoms,  H2 (v)  molecules,  F2  molecules,  HF(v'),  where  v  and 
v'  indicate  the  vibrational  quantum  states  of  H2  and  HF, 
respectively.  The  model  considers  the  ranges  of  values 
0  £  v  <_  2  and  0  v"  8,  so  that  a  total  of  16  chemical 
species  are  considered.  '  The  system  is  assumed  to  be  spatially 
homogeneous.  Over  time-scales  of  interest  it  operates  adi- 
abatically  and  at  constant  volume.  Because  of  the  chemical 


energy  release  the  temperature  is  not  constant.  Except  for 
the  Ar  atoms,  whose  number  is  constant  in  time,  the  other 
15  chemical  species  have  concentrations  which  vary  in  time. 
Including  the  temperature  variation,  the  system  is  described 
by  16  coupled  time-dependent  equations,  the  temperature  equa¬ 
tion  being  derived  from  considerations  of  energy  conservation. 
For  the  temperature  and  for  the  masses  (or  number  of  moles), 
of  all  species  except  HF(v“)  the  time  dependent  equation  is 
an  ordinary  differential  equation.  The  equations  for  the 
HF(v")  molecules  are  of  one  of  two  alternate  forms.  If  lasing 
is  not  taking  place  on  any  transition  involving  the  vibra¬ 
tional  level  v' ,  then  the  equation  for  that  level  is  an 
ordinary  differential  equation,  namely  an  equation  of  chemical 
evolution.  If  lasing  is  occurring  which  involves  the  vibra¬ 
tional  level  v"  and  an  adjacent  level  v'  +  1,  then  an  alge¬ 
braic  relation  (the  so-called  "gain  equation")  replaces  one 
of  the  chemical  differential  equations  for  the  populations  of 
state  v**  and  state  v"  +  1.  Thus  the  system  is  described  in 
terms  of  a  set  of  equations  in  which  the  equations  themselves 
change  form  in  time.  It  is  not  known  a  priori  which  equations 
are  applicable  at  a  given  time,  and  auxiliary  tests  are  re¬ 
quired  in  order  to  determine  when  changes  occur  in  the  equa¬ 
tion  system. 

An  important  parameter  of  the  system  is  the  "threshold 
gain", whose  value  depends  upon  the  reflectivity  of  the  laser 
mirrors  and  the  spacing  between  these  mirrors.  The  larger 


the  value  of  this  threshold  gain,  the  less  the  system  is  able 
to  lase.  For  sufficiently  high  threshold  gain,  lasing  is  com¬ 
pletely  suppressed.  In  this  case,  the  equations  describing 
the  system  become  completely  "chemical”  in  form,  and  consist 
entirely  of  coupled  ordinary  differential  equations.  The 
boundary  condition  on  these  equations  are  fixed  by  the  initial 
chemical  composition  of  the  system,  the  initial  pressure  and 
temperature,  and  (importantly)  by  the  number  of  fluorine  atoms 
produced  by  the  initiating  flashlamp  discharge. 

The  case  of  complete  suppression  of  lasing  action, 
which  is  usually  termed  the  “zero-power"  case,  is  therefore 
relatively  simple  to  treat,  and  it  is  also  of  considerable 
interest.  In  particular,  it  is  possible  to  study  the  timewise 
variation  of  the  populations  of  the  HF(v")  levels.  From  this 
information  one  can  calculate  the  gains  between  adjacent 
levels  and  determine  the  times  at  which  lasing  would  have 
been  initiated  had  the  threshold  gain  been  adjusted  to  some 
specific  value. 

The  model  is  characterized  by  a  very  large  number  of 
parameters,  virtually  all  of  which  are  known  with  poor  pre¬ 
cision.  The  most  significant  of  these  are  the  68  rate  co- 

* 

efficients,  one  per  reaction,  and  the  initial  conditions  of 

the  systems.  Preliminary  study  of  the  system  suggests  that 

For  each  reaction  there  are  two  rate  coefficients,  one  for¬ 
ward  and  one  reverse.  However,  appeal  to  equilibrium  con¬ 
siderations  shows  that  only  one  of  the  two  is  independent, 
and  that  the  uncertainty  in  one  stands  in  a  fixed  relation  to 
the  uncertainty  in  the  other.  Thus,  we  can  speak  of  there 
being  one  (independent)  rate  coefficient  per  reaction. 


energy  release  the  temperature  is  not  constant.  Except  for 
the  Ar  atoms,  whose  number  is  constant  in  time,  the  other 
15  chemical  species  have  concentrations  which  vary  in  time. 
Including  the  temperature  variation,  the  system  is  described 
by  16  coupled  time-dependent  equations,  the  temperature  equa¬ 
tion  being  derived  from  considerations  of  energy  conservation. 
For  the  temperature  and  for  the  masses  (or  number  of  moles) 
of  all  species  except  HF(v")  the  time  dependent  equation  is 
an  ordinary  differential  equation.  The  equations  for  the 
HF(v")  molecules  are  of  one  of  two  alternate  forms.  If  lasing 
is  not  taking  place  on  any  transition  involving  the  vibra¬ 
tional  level  v" ,  then  the  equation  for  that  level  is  an 
ordinary  differential  equation,  namely  an  equation  of  chemical 
evolution.  If  lasing  is  occurring  which  involves  the  vibra¬ 
tional  level  v"  and  an  adjacent  level  v"  +  1,  then  an  alge¬ 
braic  relation  (the  so-called  "gain  equation")  replaces  one 
of  the  chemical  differential  equations  for  the  populations  of 
state  v"  and  state  v"  +  1.  Thus  the  system  is  described  in 
terms  of  a  set  of  equations  in  which  the  equations  themselves 
change  form  in  time.  It  is  not  known  a  priori  which  equations 
are  applicable  at  a  given  time,  and  auxiliary  tests  are  re¬ 
quired  in  order  to  determine  when  changes  occur  in  the  equa¬ 
tion  system. 

An  important  parameter  of  the  system  is  the  "threshold 
gain",  whose  value  depends  upon  the  reflectivity  of  the  laser 
mirrors  and  the  spacing  between  these  mirrors.  The  larger 
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the  value  of  this  threshold  gain,  the  less  the  system  is  able 
to  lase.  For  sufficiently  high  threshold  gain,  lasing  is  com¬ 
pletely  suppressed.  In  this  case,  the  equations  describing 
the  system  become  completely  "chemical"  in  form,  and  consist 
entirely  of  coupled  ordinary  differential  equations.  The 
boundary  condition  on  these  equations  are  fixed  by  the  initial 
chemical  composition  of  the  system,  the  initial  pressure  and 
temperature,  and  (importantly)  by  the  number  of  fluorine  atoms 
produced  by  the  initiating  flashlamp  discharge. 

The  case  of  complete  suppression  of  lasing  action, 
which  is  usually  termed  the  "zero-power"  case,  is  therefore 
relatively  simple  to  treat,  and  it  is  also  of  considerable 
interest.  In  particular,  it  is  possible  to  study  the  timewise 
variation  of  the  populations  of  the  HF(v'’)  levels.  From  this 
information  one  can  calculate  the  gains  between  adjacent 
levels  and  determine  the  times  at  which  lasing  would  have 
been  initiated  had  the  threshold  gain  been  adjusted  to  some 
specific  value. 

The  model  is  characterized  by  a  very  large  number  of 
parameters,  virtually  all  of  which  are  known  with  poor  pre¬ 
cision.  The  most  significant  of  these  are  the  68  rate  co~ 

* 

efficients,  one  per  reaction,  and  the  initial  conditions  of 

the  systems.  Preliminary  study  of  the  system  suggests  that 

For  each  reaction  there  are  two  rate  coefficients,  one  for¬ 
ward  and  one  reverse.  However,  appeal  to  equilibrium  con¬ 
siderations  shows  that  only  one  of  the  two  is  independent, 
and  that  the  uncertainty  in  one  stands  in  a  fixed  relation  to 
the  uncertainty  in  the  other.  Thus,  we  can  speak  of  there 
being  one  (independent)  rate  coefficient  per  reaction. 


1 


only  13  of  the  68  rate  constants  need  be  studied  by  sensi- 
fJ  tivity  analysis,  and  that  of  the  initial  conditions,  only 

the  initial  concentrations  of  F  atoms  need  be-  studied.  We 
therefore  limited  the  sensitivity  analysis  to  14  parameters 
of  the  system. 

The  13  rate  constants  which  we  used  in  our  sensitivity  analysis  are 

detai'C-d  in  Table  4.7.  The  other  rate  constants  were  assumed  to  stand 

in  certain  fixed  relations  to  those  listed.  The  basis  for 

[13] 

this  premise  has  been  discussed  by  Cohen. 


Jhe  rate  coefficients  listed  in  . Table  4._7  are  temperature  - 


dependent. 


Broadly  speaking,  the 


values  in  Table  4.7  are  more  reliable  near  300°K  (where  they 
have  been  studied  experimentally)  than  at  higher  temperatures, 
where  numerical  evaluation  is  based  upon  extrapolation  of 
measurements  beyond  the  range  of  experimental  study.  In 
principle,  sensitivity  analysis  could  study  the  separate  in¬ 
fluences  of  uncertainties  in  the  absolute  value  of  the  rate 
coefficients  at  a  fixed  temperature  within  the  experimental 
range  (i.e.,  a  pre-exponential  or  temperature  independent 
factor) ,  and  the  influence  of  uncertainties  in  the  temperature 
dependent  part  of  the  rate  constant  (i.e.,  the  activation 
energy) .  Such  a  separation  into  temperature  independent  and 
temperature  dependent  uncertainties  was  not  undertaken  in  the 
referenced  study.  Instead,  at  the  suggestion  of  N.  Cohen^^ 
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SOME  RATE  COEFFICIENTS  OF  THE  HF  CHEMICAL  LASER^*^] 


Index  Reaction  Rate  Coefficient* 


1 

H  +  H  +  Mj  t  H2  (0)  + 

1018  T'1*0 

2 

F,  +  M,  :  F  +  F  +  M, 

2  4  4 

5.0  x 

1013  exp  (-17765 A) 

3 

HF  (v)  +  M,  J  H  +  F  +  M. 

4  4 

1.2  x 

1019  t"1,0  axp  (-68334A) 

7 

F  +  H2(0)  Z  HF(0)  +  H 

9.0  x 

1012  exp  (-805A) 

11 

HF(4)  +  H  Z  H2  (0)  +  F 

1.0  x 

1012  T0-67 

14 

H  +  F2  Z  HF(0)  +  F 

6.0  x 

1012  exp  (-1208A) 

21 

H2(l)  +  Mj^  Z  1^(0)  +  Mj^ 

2.5  x 

10*4  t4.3 

23 

HF(1)  +  Mg  Z  HF(0)  +  Mg 

9.0  x 

108  T1*3 

31 

HF(1)  +  Mg  Z  HF(0)  +  Mg 

5.0  x 

107  T1,3  +  1.0  x  1016  t"1*43 

39 

HF(1)  +  Mg  ♦  HF  (0)  +  Mg 

1.3  x 

10“2  T3,6 

47 

2HF  (1)  Z  HF  (0)  +  HF  (2) 

4.0  x 

105  T2*2 

54 

HF  (1)  +  HF(2)  Z  HF  (0)  +  HF(3) 

1.3  x 

103  T2*8 

60 

HF (1)  +  HF (3)  Z  HF (0)  +  HF(4) 

6.0  x 

10~2  t3.9 

* 

Nominal  values;  Units:  T  in  °K,  time 

in  sec, 

volume  in  an^,  mass  in  moles 

the  M^'s  denote  catalytic  species 

Mj^  =  H,  F,  Ar,  HF(0),  HF(8) 

Mj  *  20*H,  F,  Ar,  HF(0) ,  HF(8),  2.5*112(0),  2.5*H2(1), 

*  2. 5*^(2),  F2 

Mj'F. 

M4  -  H,  F,  Ar,  HF (0) ,  HF (8) ,  H2(0),  Hjd),  1^(2),  ?2 
Mg  -  H,  Ar,  Hj (0) ,  H^ (1) ,  H2 (2) ,  F2 
Mg  *  HF(0) ,  HF(8) 

Coefficients  such  as  "20*"  In  the  list  of  catalytic  species  denote  that  the 
species  concentration  is  weighted  by  this  factor  in  computing  the  reaction  r 

-6^- 
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the  rate  coefficients  were  assumed  to  be  undertain  to  within  a 
multiplicative  factor  or  5,  independent  of  temperature.  \ 

At  every  temperature  each  rate  coefficient  k-j  was  assumed 

to  have  a  uniform  probability  distribution  of  its  logarithm 

within  bounds  given  by 

log  -  log  5  <  log  ki  <  log  k±  +  log  5.  (4.13) 

where  k^  is  the  nominal  value  of  k^,  i.e.,  the  value  tabulated 
in  Table  4.7.  We  note  in  passing  that  for  such  a  distribution, 
the  transformation  function  G^  of  Eq.  (3.1)  is  given  by 

Gi(x)  =  |  Sin-1  (x)  ,  (4.14) 

where  Sin-1  is  the  principal  value  of  the  inverse  sine  function. 

The  number  of  fluorine  atoms  produced  at  time  t  *  0  by 

flashlamp  discharge  is  similarly  uncertain,  by  about  a  factor 

ns] 

of  two  »i.e.,  for  this  parameter  log  5  is  replaced  by 
log  2  in  Eq.  (4.13)  . 

The  differential  equation  system  which  describes  the 
laser  at  zero  power  was  programmed  and  run.  Initial  condi¬ 
tions  of  temperature  and  composition  were  as  indicated  in 
Table  4'.  8. 

The  laser  model  was  integrated  from  zero  time  out  to 

a  real  time  of  4.0  microseconds.  With  fourteen  parameters 

* 

it  was  necessary  to  carry  out  907  such  integrations.  Data 

This  is  based  on  the  rule  N  =  2  w  -  7  mentioned  previously, 
In  Section  III  and  not  on  the  rule  N  *  2w  X  +  1 . 


TABLE  4 .8 


INITIAL  CONDITIONS  FOR  LASER  STUDY 


Temperature  300°K 
Species  Concentrations 

Ar  (inert  diluent) 

h2(0) 

H2(l) 

H2(2) 

H 

HF  (v*0)  (v"  *  0  to  8) 

F  +  2F2 

F  (nominal  value) * 

** 

F2  (nominal  value) 


4.704 

x  10“5  mole  era”3 

9.407 

X 

10'7  mole  cm”3 

1.119 

X 

10”15  mole  cm”3 

4.321 

X 

-24  -3 

10  mole  cm 

0.0 

mole  cm”3 

0.0 

mole  cm”3 

1.975 

X 

10”®  mole  cm”3 

4.704 

X 

—8  —3 

10  mole  cm 

1.928 

X 

10"®  mole  cm”3 

*  Actual  initial  value  ranges  from  -50%  to  100%  of 
tabulated  value  (see  text) . 

**  Actual  initial  value  dependent  upon  initial  F  atom 
concentration,  such  that  sum  F  +  2F2  is  tabulated 
value. 


on  populations  of  vibrational  states,  chemical  concentrations, 
temperature,  and  other  variables  was  stored  at  each  100  nano¬ 
second  real  time  interval.  -From  the  populations  N(v)  of 
adjacent  vibrational  levels  of  HF,  it  is  simple  to  calculate 
the  gains.  The  gain  a(v,J)  for  the  transition  HF(v+l,  J-l) 

-*•  HF  (v,  J)  (J  is  the  rotational  quantum  number  of  the  lower 
vibrational  state)  is  expressed  by 

c*(v,J)  -  ^  wc(v,J)$  (v,J)B(v,J)  j|jrr  N(v+1,J-1)-N(v,J)J 

(4.15) 

where  h  is  Planck's  constant,  is  Avogadro's  number,  u>c 

is  the  wavenumber  of  the  transition,  $  (v,  J)  is  the  line  pro- 

* 

file  at  line  center,  B  (v,J)  is  the  Einstein  coefficient  for 
absorption,  and  N{v,J)  is  the  population  of  the  v,J  rotation- 
vibration  level  of  HF.  The  assumption  of  thermal  equilibrium 
of  the  rotational  states  implies  that  N(v,J)  satisfies 

N (v, J)  -  N(v)  (2J+l)exp(-hCEj/kT)/Q^(T)  (A16) 

where  k  is  Boltzmann's  constant,  T  is  the  absolute  temperature 
(assumed  the  same  for  rotation  and  translation) ,  Ej  is  the  „ 
rotational  energy  of  the  v,J  level  relative  to  the  v,0  level, 
and  Q^(T)  is  the  rotational  partition  function; 

Qr(T)  "  S  (2J+l)exp(-hcEj/kT).  (‘4.17) 

J  *  -  • 

In  evaluating  the  gain  as  a  function  of  time  it  is  necessary 
to  search  over  the  rotational  quantum  number  J  to  find  that 
The  Doppler  profile  was  assumed  in  the  model. 
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J  value  which  maximizes  a(v,j)  as  a  function  of  time.  The 
reason  for  this  is  that  only  one  such  level  can  lase  at  any 
given  time  within  each  vibrational  band,  and  it  is  the  level 
for  which  ot(v,J)  is  maximum  that  actually  lases.  This  J 
value  can  shift  in  time.  The  search  procedure  is  simple  and 
straightforward . 

We  will  describe  the  results  of  the  sensitivity 

analysis  for  the  v  =  2-*-l,  v  =  3-*-2,  and  v  =  4-**3  bands.  The 

r  9] 

report  by  Levine1-  '*  analyzed  many  other  variables  in  addition 
to  these  three,  and  the  reader  can  refer  to  this  report  for 
additional  discussion.  The  three  variables  we  will  describe 
here  suffice  to  illustrate  the  technique  of  application  of 
sensitivity  analysis. 

Figure  ~4~  shows  the  zero-power  gain  as  a  function  of 
time  for  the  v  =  2-*l,  v  *  3-*2,  and  v‘  =  4-*-3  transitions.  Both 
the  nominal  values  and  the  mean  values  (averaged  over  the  dis¬ 
tribution  of  all  parameters)  are  shown.  Except  at  times 
less  than  1  usee  following  initiation  the  mean  and  nominal 
values  differ  considerably,  both  for  the  v  =  2-*-l  and  v  *  3-*2 
transition.  For  the  v  =  4-*3  transition,  the  mean  and  nominal 
values  are  not  in  good  agreement  even  at  times  as  short  as 
0.2  usee.  From  this  we  can  infer  that  the  variance  of  these 
transitions  will  be  large,  which  is  confirmed  in  Figure  5  , 

where  we  plot  the  coefficients  of  variation.  From  these  two 
figures  we  see  that  the  predictive  ability  of  the  model  is 

poor  as  a  consequence  of  the  parameter  uncertainties. 


The  question  as  to  which  parameter  uncertainties  cause  this 
high  variance  is  answered  graphically  in  Figures  6,  7, 
and  8.  These  three  figures  display  the  partial  variances 
for  the  v  =  2-*l,  v  =  3-*2,  and  v  =  4-*3  gains,  respectively,* 
Reference  to  Figure  6  shows  that  for  the  v  *  2-*-l 
transition  the  variance  at  early  times  is  due  mainly  to  un¬ 
certainties  in  the  initial  F  atom  concentration  [F] ^  and  the 
rate  coefficient  for  process  F  +  H2(0)  Z  HF(O)  +  H.  Since 
the  coefficients  for  the  process  F  +  H2 (0)  Z  HF(v) ;  +  H 
(v  =  1,  2,  3)  are  proportional  to  this  rate  (c.f..  References 
9/and 14)  this  is  equivalent  to  stating  that  calculation  of  the 
gain  at  early  times  is  sensitive  to  the  rates  at  which  re¬ 
action  of  F  with  H2 (0)  populates  the  excited  levels,  as  we 
should  anticipate. 

At  longer  times  the  sensitivity  coefficient  for  the 

collisional  deactivation  process  HF(v)  +  HF(v')  Z  HF(v-l)  +  HF(v') 

becomes  the  dominant  cause  for  uncertainty  in  the  computed 

values  of  the  gains.  At  still  longer  times  the  uncertainties 

in  the  rate  coefficient  for  the  process  H  +  F2  Z  HF(0)  +  F, 

and  in  the  initial  concentration  of  F  atoms,  become  the  dominant 

sources  of  uncertainty.  The  reason  for  the  importance  of  the 

latter  parameter  at  late  times  is  related  to  the  fact  that  the 

T  91 

laser  at  zero-power  operates  adiabatically.  The  reason  for 

the  importance  of  the  foimer  parameter  at  late  times  is  that 

the  collection  of  reactions  H  +  F2  Z  HF(v)  +  F  tends  to 
_  — — 

The  curves  only  show  those  partial  variances  which  are  large. 
Parameters  whose  partial  variances  are  negligible  are  not  shown, 
so  as  to  avoid  cluttering  the  figures. 
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repopulate  the  excited  levels  HF(v)  (v  1)  following  their 
initial  depopulation.  T  he  rate  coefficients  for  this  process 
with  v  >_  1  are  taken  to  be  proportional  to  the  rate  coefficient 
for  this  process  with  v  =  0,  c.f..  References  9  andl3‘). 

Reference  to  Figure  7„.  shows  that  the  variance  in  the 
computed  zero-power  gain  at  the  v  =  3-*- 2  transition  is  due  to 
substantially  the  same  reasons  as  in  the  case  v  =  2-*l,  with 
small  additional  influences.  The  latter  correspond  to  the 
processes  HF(v)  +h;f+  H2(0)  (v  =  4,  5,  6)  and  2HF(v)  * 
HF(v-l)  +  HF (v+1) .  The  former  process,  while  it  does  not 
directly  involve  the  states  HF(2)  and  HF(3),  nonetheless  in¬ 
fluences  these  states  indirectly,  since,  e.g.,  removal  of 
molecules  from  HF(4)  by  collision  with  H  eliminates  molecules 
in  HF (4)  as  a  source  of  molecules  in  HF(3)  via  the  process 
HF  (4 )  +  HF  (v)  t  KF(3)  +  HF  (v)  . 

Reference  to  Figure  8^  shows  a  feature  unique  to  the 
zero-power  gain  on  the  v  =  4-*3  transition.  The  process 
HF (v)  +  H  t  H2(0)  +  F  (v  =  4,  5,  6)  is  the  dominant  cause  of 
variance  in  the  computed  gain  at  times  beyond  1.0  usee. 
Furthermore,  as  can  be  seen  in  Figure  5,  j  the  variance  for 
the  gain  on  this  transition  is  larger  than  that  for  the 
other  transitions.  It  is  striking  also  to  note  that  the 
computed  gain  on  this  transition  is  very  small  (see  Figure 
4)j  a  fact  which  can  be  related  to  a  very  low  population 
for  the  v  *  4  level. It  thus  becomes  possible  to  sug¬ 
gest,  as  a  consequence  of  the  sensitivity  analysis,  that  the 


anomalously  small  computed  gain  cn  the  v  =  4-*  3  transition 

may  be  due  to  too  large  a  nominal  value  for  the  rate  coefficient  of  the  process 

HF(4)  +  H  J  1^(0)  +  F*  For  ^urt^*er  discussion  of  this 

19  i 

point/  we  refer  the  reader  to  the  report  by  Levine. 1  1 

In  concluding  this  example,  it  is  worth  noting  that 
the  results  of  the  sensitivity  analysis  show  that  the  pre¬ 
dictions  of  the  model  are  sensitive  only  to  five  of  the 
14  parameters  studied.  Thus  we  find  that  the  large  un¬ 
certainties  in  the  remaining  nine  parameters  do  not  influence 
the  predictions.  In  view  of  this,  we  can  conclude  that  for 
all  practical  purposes  those  processes  for  which  the  sensi¬ 
tivity  coefficient  is  small  can  be  omitted  from  the  model 
entirely.  That  is  to  say,  a  simplified  five  reaction  model 
would  suffice  to  describe  the  HF  laser  at  zero  power. 
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A  SOCIOECONOMIC  MODEL  (WORLD  II) 


Forrester  has  written  extensively  on  the  use  of  mathe¬ 
matical  models  as  descriptions  of  industrial  and  social  pro¬ 
cesses.^*^  World  II  an  example  of  a  model  of  this 

type.  It  consists  of  five  coupled  differential  (or  dif- 

* 

ference)  equations  describing  the  time  dependence  of  five 
variables:  population,  natural  resources,  capital  invest¬ 
ment,  pollution,  and  the  fraction  of  capital  which  is  in 
agriculture.  The  model  utilizes  a  large  number  of  parametric 
relations  to  relate  the  time  rates  of  change  of  the  dependent 
variables  to  current  values  of  these  variables.  Many  of 
these  relations  consist  of  tabulations,  so  that  each  table 
entry  constitutes  a  parameter.  As  such,  the  model  contains 
several  hundred  parameters,  each  of  which  contains  to  some 
extent  biases  of  the  model's  author.  In  any  event,  none  of 
these  parameters  are  precisely  known,  and  these  parametric 


* 

Typically  the  models  developed  by  Forrester  consist  of 
coupled  difference  equations  in  a  set  of  dependent  variables 
with  time  as  independent  variable.  The  change  xn  -  xn„i  in 
a  typical  dependent  variable  over  one  time  interval 
At  *  nAt  -  (n-l)At  is  usually  expressed  as  a  forward  dif¬ 
ference,  i.e.,  x^  -  xn_i  is  assumed  to  depend  upon  the  set 
of  dependent  variables  txn_^},  plus  in  some  case  independent 
driving  functions,  all  evaluated  at  time  t  *  (n-l)At.  In 
the  limit  At  0  the  system  of  equations  would  reduce  to  a 
system  of  ordinary  differential  equations.  As  compared  to 
the  differential  equations  so  derived,  Forrester's  difference 
equations  can  be  viewed  as  being  the  forward  Euler  difference 
approximation  to  the  differential  equation.  As  such,  the 
possibility  arises  that  the  difference  equations  may  be 
relatively  unstable  compared  to  the  differential  equations, 
and  that,  therefore,  the  difference  form  may  generate  solu¬ 
tions  which  depend  sensitively  upon  the  time  step  employed. 

To  our  knowledge,  no  specific  stability  analysis  of  Forrester's 
World  II  model  has  been  performed,  but  the  numerical  solu¬ 
tions  show  no  evidence  of  artificial  mathematical  instabilities 
induced  by  the  difference  scheme  used. 


uncertainties  will  influence  the  prediction  of  the  model.  We 
selected  a  small  subset,  nine  in  number, of  the  parameters 
and  examined  the  sensitivity  of  the  model  to  assigned  uncer¬ 
tainties  in  these  parameters.  The  parameters  we  studied 
are  listed  in  Table  4.9,  along  with  the  "identifiers"  used 
in  the  original  study. 

The  first  four  parameters  listed  in  Table  4.9  are 
annual  rates  of  growth  or  decay.  The  last  of  these  para¬ 
meters,  CIDN,  basically  is  identical  to  the  rate  of  deprecia¬ 
tion  of  capital.  The  fifth  parameter,  "normal  pollution",  is 
a  multiplier  which  gives  the  rate  of  generation  of  pollution 
per  unit  of  population,  the  unit  of  pollution  being  given 
by  definition  in  terms  of  a  standard  unit  defined  so  that  the 
per  capita  global  pollution  in  1970  is  one  unit. 

The  sixth  parameter,  BETA,  is  related  to  the  rate  of 
absorption  of  pollution,  i.e.,  its  conversion  to  harmless 
form,  by  the  global  ecosystem.  Forrester  assumes  aquasilinear 
relation  between  absorption  time  and  level  of  pollution  (see 
Figure  3-15,  p.  57,  of  Reference  17; .  For  such  a  case,  the 
rate  of  change  of  pollution  p,  would  be  given  by 

dp/dt  =  -p/a(p)  +  source  terms  (4.19) 

where  the  characteristic  decay  time  a  is  itself  a  function  of 
p.  We  have  scaled  Forrester's  estimate  of  a(p)  by  a  constant 
factor  S ,  i.e.,  we  use 

dp/dt  =  -p/8a(p)  +  source  terms. 


(4.20) 


TABLE  4  > 9 


PARAMETERS  OF  WORLD  II  MODEL  STUDIED  BY  SENSITIVITY  ANALYSIS 


Identities 

Nominal  Value 

1. 

BFN 

Normal  birth  rate 

0.04  year”1 

2. 

DHN 

Normal  death  rate 

0.028  year”1 

3. 

CIGN 

Normal  capital  investment 
growth  rate 

0.05  year”1 

4. 

CIDN 

Normal  capital  investment 
decay  rate 

0.025  year”1 

5. 

POLN 

Normal  pollution 

1.0 

6. 

BETA* 

Scale  factor  for  pollution 
absorption  time 

1.0 

7. 

CIAFT 

Capital  fraction  in  agriculture 

15.0  years 

8. 

C* 

Coefficient  of  CFIFR  in  rate  of 
change  of  CIAF  (see  text) 

1.0 

9. 

NRI 

Initial  inventory  of  natural 
resources 

9.0  x  10U 

* 

Not  defined  in 

original  world  II  model.  See  text 

for  explanation. 
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and  our  sixth  parameter  is  this  scale  factor  8  (BETA) . 

Values  of  BETA  in  excess  of  unity  correspond  to  ecosystem 
recovery  slower  than  that  assumed  by  Forrester;  values  of 
BETA  less  than  unity  correspond  to  ecosystem  recovery 
faster  than  that  assumed  by  Forrester. 

The  seventh  parameter  is  a  relaxation  time  which  fixes 
the  rate  at  which  capital  investment  can  be  shifted  into  or 
out  of  agriculture  in  response  to  the  variations  in  demand 
for  agricultural  capital  (e.g.,  tractors,  harvesting  equip¬ 
ment)  ,  this  demand  also  being  dependent  upon  the  instantaneous 
supply-demand  relation  for  food. 

The  eighth  parameter  is  similar  in  character  to  the 
sixth.  Forrester  assumes  that  the  rate  of  change  of  the 
fraction  of  capital  investment  devoted  to  agriculture  is  re¬ 
lated  to  the  per  capita  food  supply  (see  Figure  3.13,  p.  59 
of  Reference  1 7 » ,  with  larger  food  availability  leading  to 
decreasing  rate  of  capital  investment  in  agriculture,  and 
vice  versa.  His  assumed  relation  is  roughly  exponential, 
with  all  new  capital  invested  in  agriculture  in  the  limit 
of  zero  food  availability;  with  30  percent  of  new  capital  de¬ 
voted  to  agriculture  when  the  per  capita  food  availability 
(the  "food  ratio")  is  at  a  "satisfactory"  level  (unity  by 
definition) /  and  '9  percent  of  new  capital  when  food  supply 
is  twice  per  capita  need.  We  can  express  this  relation  as 
f  *  exp  (-1.204  r)  where  f  is  the  capital  fraction  (CIAF) 
indicated  by  the  food  ratio  (CFIFR)  and  r  is  the  food  ratio. 
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In  our  analysis,  we  generalize  this  relation  to  read 
f  =  Cexp(-1.204r) ,  so  that  values  of  C  in  excess  of  unity 
correspond  to  increasing  the  rate  at  which  new  capital  is 
switched  into  agriculture  in  response  to  food  shortages, 
as  compared  to  the  standard  rate  assumed  by  Forrester,  and 

values  of  C  less  than  unity  correspond  to  decreasing  this 

* 

rate  as  compared  to  Forrester's  assumption. 

The  last  parameter  is  the  inventory  of  natural  re¬ 
sources  in  the  initial  year  of  the  model,  in  arbitrary  re¬ 
source  units.  The  nominal  value  of  this  inventory  is  what 
is  presumed  available  in  the  year  1900,  with  a  per  capita 

need  of  1  unit/year.  This  supply  then  would  suffice  for 

g 

250  years  for  a  (fixed)  population  of  3.6  *  10  individuals. 

A  major  criticism  of  World  II  has  been  that  its  re¬ 
source  inventory  is  finite  and  irreplacable.  In  such  a 
case  this  initial  value  must  become  the  dominant  parameter 
at  sufficiently  long  times.  Rather  than  try  to  introduce 
renewable  resources,  technological  innovation,  or  other 
mitigating  concepts,  we  have  chosen  to  analyze  World  II  as 
is.  Therefore,  we  will  inevitably  be  confronted  with  a 
world  in  which  everyone  dies  in  the  long  term.  At  issue, 
however,  is  the  question  of  the  reliability  of  World  II  on 
a  relative  basis  in  the  intervening  time. 

* - 

At  first  sight  it  might  appear  that  when  C  >  1  and  the 
shortage  of  food  is  acute,  the  system  will  "try"  to  devote 
more  than  100  percent  of  all  new  capital  to  agriculture. 

Needless  to  say,  it  cannot  succeed  in  doing  so,  as  careful 
analysis  shows.  But  it  does  mean  that 

net  rate  of  response  to  "starvation"  is  better  than  in  the 
standard  case. 

-TV 


One  auxiliary  variable  is  important.  Forrester  de¬ 
fines  a  "quality-of-life"  function  as  the  product  of  four 
factors,  each  dependent  upon  one  of  the  dependent  variables: 
population,  food,  capital  investment  per  capita  (termed 

"material"),  and  pollution.  There  are  several  arbitrary 

* 

aspects  to  his  definition,  such  that  the  numerical  value 
is  best  viewed  only  as  qualitative  and  not  quantitative. 
Nonetheless,  to  lay  stress  on  the  difficulty  of  predicting 
qualify  of  life,  we  have  treated  this  variable  as  quantita¬ 
tive,  and  have  subjected  its  predicted  values  to  sensitivity 
analysis. 

** 

We  obtained  a  computer  program  for  World  II,  and 
applied  our  method.  Except  for  the  nine  parameters  listed 
in  Table  4.9,  all  parameters  were  given  the  values  used  by 
Forrester  for  his  "base  case".  For  the  nine  parameters  of 
Table  4.9,  the  nominal  values,  as  given  in  the  table,  cor¬ 
respond  to  Forrester's  base  case.  We  attributed  to  each  of 

*  We  are  tempted  to  suggest  that  a  Delphi  procedure  might  pro¬ 
vide  the  most  rational  way  of  determining  a  consensus  defini 
tion  of  quality-of-life. 

** 

We  are  indebted  to  Professor  W.  E.  Schiesser,  Department  of 
Chemical  Engineering,  Lehigh  University,  for  supplying  us 
with  this  program.  This  program  uses  second  order  forward 
Runge-Kutta  integration  instead  of  first  order  forward 
Euler  differencing.  It  is  possible  to  prove  that  this 
Runge-Kutta  method  has  wider  stability  bounds  than  the  for¬ 
ward  Euler  method,  so  that  with  equal  or  smaller  time-steps 
than  those  used  by  Forrester,  we  are  assured  that  our  solu¬ 
tions  display  the  same  stability  character  as  Forrester's 
original  equations.  We  in  fact  used  smaller  time- 
steps,  so  that  numerically  we  always  generated  solutions 
with  the  same  stability  characteristics  as  in  the  original 
work. 
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these  nine  parameters  an  uncertainty  of  +  20  percent#  with 
a  uniform  distribution  of  probability  within  this  interval. 

In  view  of  the  difficulty  in  precisely  ascertaining  values 
for  these  parameters,  it  appears  to  us  that  attribution  of  a 
+20  percent  error  understates  the  real  situation. 

The  output  (dependent)  variables  we  analyzed  were 
population#  pollution#  capital  investment#  natural  resources# 
and  quality  of  life.  In  this  paper  we  only  discuss  the 
first  and  last  of  these. 

Figures  9~  and  _10.  show  our  results  for  population. 

In  Figure  9  J  we  see  that  the  nominal  and  mean  values  of 
population  are  virtually  coincident  until  about  Year  1976 

(they  also  coincide  fairly  well  with  real  data)  but  that 

* 

beginning  about  Year  1976  a  divergence  sets  in.  This  diver¬ 
gence  reaches  an  (absolute)  maximum  about  Year  2040#  at  which 

g 

time  the  mean  value  is  about  one  billion  (1  *  10  )  persons 
lower  than  the  nominal  value.  The  corresponding  value  of 
the  coefficient  of  variation  also  is  shown;  it  reaches  a 
maximum  value  of  about  0.4  in  Year  2040. 

The  sensitivity  coefficients  for  population  are  shown 
in  Figure  10.  i  The  graph  exhibits  only  four  of  the  nine  para¬ 
meters  because  for  the  other  five  the  coefficients  are  very 

*We  see  here  a  remarkable  example  of  the  well-known  fact 
that  models  do  a  better  'job  of  confirming  the  past  than  of 
predicting  the  future! 


small,  and  effectively  zero.  It  is  interesting  to  note  that 
in  early  years,  i.e.,  pre  1930,  the  birth  and  death  rates, 
which  are  direct  influences  of  population,  are  the  dominant 
parametric  influence  in  population.  However,  at  later  times, 
the  normal  capital  investment  growth  rate  emerges  as  the  most 
important  parameter,  and  the  variance  in  population  prediction 
is  dominated  by  the  uncertainty  in  this  investment  rate. 

The  peculiar  effect,  near  Year  2010,  where  nearly  all 
of  the  sensitivity  coefficients  fall  to  zero  arises  from  an 
inadequacy  in  our  procedure  as  applied  to  this  one  example; 
at  the  time  we  did  this  calculation  (mid  1974)  we  were  not 
yet  aware  of  the  significance  of  the  harmonics  in  the  defini¬ 
tion  of  the  sensitivity  coefficients,  so  that  we  were  com¬ 
puting  coefficients  by  calculating  the  fundamental  terms  of 
Eq.  (2.25)  only.  In  such  a  calculation,  it  is  possible  to 
"lose"  some  of  the  variance  when  the  coefficient  of  the 
Fourier  fundamental  undergoes  a  change  in  sign  (c.f..  Section  II 
1  for  a  discussion  of  this  possibility).  This  circumstance 
applies  for  population  in  Year  2010.  It  would  be  interesting 
to  repeat  our  calculation  using  the  fuller  definition,  Eq. 

(2.25) ,but  we  have  not  had  the  opportunity  to  do  so.  None¬ 
theless,  it  seems  clear  to  us  that  the  uncertainty  in  the 
normal  capital  investment  growth  rate  is  the  key  parameter 
which  leads  to  variance  in  the  predicted  population. 

Figures  11  'and  12  show  our  analysis  of  the  auxiliary 
variable  quality  of  life.  In  Figure  12  "  we  see  that  the 


mean  and  nominal  values  of  this  variable  are  virtually  co¬ 
incident  until  about  Year  2005,  but  that  following  this  year 
the  mean  value  rises  to  a  dramatic  new  level  whereas  the 
nominal  value  continues  to  decay.  After  about  Year  2045 
the  mean  value  rapidly  decays.  The  coefficient  of  variation 
is  small  until  Year  2005,  but  then  rises  dramatically  over 
the  next  15  years,  following  which  it  too  decays.  The  long 
term  decay  (i.e.,  after  Year  2050)  of  the  mean  and  nominal - 
values  and  the  coefficient  of  variation  all  are  traceable 
to  the  finite  resource  inventory  in  the  model.  The  influence 
of  the  latter  only  is  beginning  to  appear  at  about  Year  2070, 
but  is  clear  enough  that  with  an  inventory  capable  of  sus- 

q 

taining  3.6  *  10  people  for  250  years,  and  populations  as 
exhibited  in  Figure  9,  depletion  of  resources  is  going  to 
be  of  overriding  significance  in  the  model  by  Year  2150  * 

1900  +  250,  i.e.,  just  off  the  figure  at  the  right. 

It  is  striking  to  see  that' the  model  is  more 
likely  to  lead  to  a  prediction  of  a  "golden  age"  in  the 
period  2005-2050  than  it  is  to  predict  a  continuing  decay, 
if  an  accounting  is  made  of  parameter  uncertainties. 

Figure  1Z  1  shows  the  sensitivity  coefficients  for  the 
quality  of  life:  birth  and  death  rate  are  the  important 
parameters  prior  to  2005;  however,  since  the  variance  is 
small  in  these  years,  it  is  not  too  significant  to  inquire 
as  to  the  causes  of  the  variance.  After  Year  2020,  the 
normal  capital  investment  growth  rate  becomes  the  most 
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important  single  parameter,  although  it  by  no  means  dominates 
quality  of  life  to  the  degree  it  dominates  population.  Still 
later  in  time,  there  is  a  transient  period  of  importance  for 
the  normal  capital  investment  decay  rate,  presumably  due  to 
the  fact  that  with  diminishing  resources  and  hence  limited 
capability  for  new  capital  formation,  it  is  the  rate  of  dis¬ 
appearance  of  capital  to  which  quality  of  life  is  most  sensi¬ 
tive. 

To  summarize  then,  the  conclusions  we  have  reached  in  our  brief 
study  of  Forrester's  model  are: 

1.  Uncertainties  as  to  the  nominal  rate  of  capital  in¬ 
vestment  is  the  dominant  parameter  which  renders  pre¬ 
dictions  of  World  II  unreliable,  at  least  over  the 
next  100  years. 

2.  To  the  extent  to  which  the  resource  inventory  is  truly 
finite,  the  initial  extent  of  this  resource  ultimately 
becomes  the  key  parameter  at  some  future  point  in  time. 

3.  The  model  is  not  capable  of  predicting  "quality  of 

life"  to  any  meaningful  degree  after  about  Year  2005. 

(18) 

More  recently,  Meadows  et  al.  1  have  developed  World  3,  an  extended 
and  enlarged  version  of  World  2,  with  many  more  variables,  parameters,  and 
equations.  They  subject  this  model  to  a  number  of  linear  sensitivity  analyses, 
on  the  basis  of  which  they  determine  the  model's  qualitative  reliability.  A 
much  more  appropriate  test  of  the  model's  reliability  would  be  via  a  nonlinear 
method  such  as  the  one  discussed  in  this  paper. 
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A  CH EH  I CAL  REACTION  MODEL 


(19) 

In  a  recent  paper'  ,  Boni  and  Penner  have  successfully  applied  our 
method  of  sensitivity  analysis  to  a  study  of  methane  oxidation  kinetics. 
Their  model  consisted  of  a  set  of  23  co-ioled  rate  equations  involving  23 
parameters,  i.e.  rate  coefficients.  These  parameters  were  varied  over  a 
+50%  range  of  uncertainties  about  their  nominal  values  kj0^.  The  output 
functions  were  the  species  concentrations  [CH^,  CH^,  CHgO,  CHO,  CO,  C02, 

H20,  0,  H,  OH]  at  different  times  [10"^  to  10-2  seconds]  after  the  initiation 
of  the  reaction.  Their  analysis  showed  "that  of  the  23  reactions,  plus  their 
inverses,  included  in  the  mechanism,  only  about  5-7  reactions  (depending  upon 
the  species  in  question)  strongly  affect  the  concentration  of  that  species". 
Their  analysis  thus  provides  an  important  example  of  how  sensitivity  analysis 
can  be  used  to  simply  complex  models  by  segregating  the  "important"  and 
"unimportant"  component  equations. 


V.  Additional  Research 


There  are  a  number  of  interesting  and  important  problems  which  have 
arisen  during  the  course  of  this  research  and  to  which  we  have  not  had  an 
opportunity  to  address  ourselves.  We  list  them  here  briefly  and  hope 
that  some  of  the  readers  of  this  review  and  future  practioners  of  this  method 
will  be  stimulated  to  carry  out  some  work  along  these  lines, 
a)  Postponement  of  Interferences: 

In  sections  2E  and  2F  above  and  in  Appendix  I,  we  have  discussed  the 
problem  of  interference  which  arises  from  the  unavoidable  use  of  a  set  of 
integer  or  rational  frequencies.  This  effect  can  be 

"postponed"  in  the  sense  that  we  can  chose  a  high  value  of  M  (see  Eq.  A1.4), 
the  order  of  interference.  As  pointed  out  in  refs.  (1)  through  (3),  how¬ 
ever,  the  larger  the  chosen  value  of  M,  the  larger  the  maximum  value  w  v 
of  the  input  frequencies  set  {w}  and  correspondingly,  the  larger  the  number 
N  of  s-space  points  required  for  the  evaluation  of  the  Fourier  amplitude. 

Thus,  a  large  value  of  M,  which  will  minimize  the  interferences,  will 
appreciably  increase  the  number  of  computations  required  for  the  calculation 
of  the  Fourier  amplitudes.  As  we  have  shown,  for  instance,  in  ref.  (3), 

806  points  in  s-space  are  required  for  the  calculation  of  the  Fourier 

amplitudes  for  a  10  parameter  system  for  M=4,  while  8,520  points  are 

* 

required  for  the  same  system  when  M=6.  Since  o>  will  of  course  increase 
with  the  number  n  of  parameters  k. ,  i=l,2,...,n,  this  problem  becomes 
particularly  serious  for  systems  with  a  large  number  of  parameters. 

One  obvious  way  to  circumvent  this  problem  completely  would  be  to  use 
a  set  of  incommensurate  frequencies  u.  Since  this  does  not  seem  a  feasible 
on  a  computer  with  a  finite  register,  the  next  best  solution  would  be  to 

*0wing  to  the  symmetry  properties  of  f's)  discussed  in  Section  (3.2),  the  number 
of  sample  points  required  are  only  one-half  of  those  listed  here. 


learn  how  to  construct  integer  frequency  sets  {w}  which  will  lead  to 


large  values  of  fl  for  reasonable  values  of  ojmax .  It  is  not  clear  how  much 

of  an  improvement  can  be  achieved  on  the  frequency  sets  already  published 
(1  2  3) 

in  our  papers'  *  ’  but  research  along  these  lines  is  clearly  desirable. 

In  this  connection,  attention  should  be  called  to  our  brief  discussion 
in  section  2F  (see  Eq.  2.20)  on  the  spacing  of  the  N  quadrature  points 
used  in  the  calculation  of  the  Fourier  coefficients.  It  may  well  be  possible 
to  increase  the  accuracy  of  the  computed  Fourier  coefficients  by  a  more 
judicious  choice  of  spacing  which  takes  account  of  the  oscillatory  properties 
of  the  output  function  f(s).  Some  work  on  this  problem  could  be  very  useful. 

b)  Expansion  of  Output  Function  f_ 

We  have  chosen  in  the  work  done  so  far  to  transform  the  output  functions, 
f  into  periodic  functions  on  (0,  2*),  see  Eq.  (2.7),  and  then  Fourier  analyse 
these  functions  to  obtain  their  Fourier  coefficients.  Interestingly,  there 
is  nothing  sacrosanct  in  expanding  f  in  terms  of  sine  and  cosine  functions. 

One  could  equally  well  expand  f(u),  the  output  function  in  u-space,  in 
terms  of  any  other  desired  or  useful  set  of  orthogonal  functions,  such  as, 
for  instance,  Hermite  polynomials.  One  would  then  have  to  establish  again 
the  connection  between  the  expansion  coefficients  and  some  "sensitivity  measure" 
as  has  been  done  above.  It  is  not  clear  whether  such  expansions  in  terms 
of  other  orthogonal  functions  would  lead  to  a  simpler  or  better  (or  worse) 
theory  of  nonlinear  sensitivity  analysis,  but  it  raises  a  question  which 
could  usefully  be  pursued. 

c)  Information  in  Harmonics  and  Combination  Frequencies 

As  discussed  in  the  previous  sections,  the  Fourier  spectrum  of  the 
output  function  contains  harmonics  and  linear  combinations  of  all  the  input 
frequencies  «  ,  i  -  1 ,2,. . . ,n.  In  our  earlier  version  of  sensitivity  analysis, 
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references  (1)  through  (3),  we  have  used  only  the  information  contained  in 

the  fundamentals  u>„  through  the  Fourier  amplitudes  B  ,  Eq.  (2.30).  In  the 
i 

partial  variance  method  outlined  in  section  b  above  we  have  made  use  also 

of  the  harmonics  of  the  fundamental  frequencies  for  the  construction  of  S 

“i 

as  in  Eq.  (2.26).  As  pointed  out  in  that  section,  Eqs.  (2.27)  through  (2.29), 

2  2 

it  is  also  possible  to  construct  higher  partial  variances  S  ,  S  , 

V“k 

etc.  corresponding  to  linear  combinations  of  the  fundamental  frequencies 
Wj,  etc.  These  higher  partial  variances  contain  increasingly  more 
detailed  information  about  the  coupling  of  sensitivity  due  to  uncertainties 
of  groups  of  parameters.  Such  information  is  clearly  of  great  importance 
in  studying  the  sensitivity  of  chemical  and  other  complex  rate  systems  since 
the  explicit  rate  laws  for  the  output  functions,  if  they  could  be  obtained, 
would  most  probably  involve  sums  and  products  of  various  parameters.  It  would 
therefore  be  most  useful  to  explore  the  construction  of  higher  partial  variances 
and  analyze  the  information  obtained  from  them  in  future  applications  of 
the  Fourier  method  of  sensitivity  analysis, 
d)  Correlated  Parameters 

In  all  the  work  described  so  far  we  have  made  the  explicit  assumption 
that  the  system  parameters  k  and  their  variation  u  over  any  desired  range 
are  uncorrelated.  By  this  we  mean  that  each  parameter  can  be  varied 
Independently  of  all  other  parameter  or,  equivalently,  that  to  each  parameter 
one  can  assign  a  range  of  uncertainty  with  its  probability  distribution 
independent  of  the  uncertainty  range  assigned  to  all  other  parameters.  This 
is  certainly  valid  for  many  physical  systems  where  the  parameters  are  indeed 
Independent  of  one  another  and  can  be  determined,  theoretically  or  experi¬ 
mentally,  independently  of  one  another. 


In  many  economic  and  social  model  systems  there  are,  however,  correla¬ 
tions  betewen  the  parameters.  Frequently  it  is  not  even  possible  to  define 
parameters  which  are  independent.  The  determinations  of  parameters  in  such 
systems  by  fitting  models  to  repeated  observations  leads  to  a  set  of 
parameters  which  are  statistically  dependent,  i.e.  the  range  of  uncertainty 
of  parameter  may  well  be  correlated  with  the  uncertainty  range  of  para¬ 
meters  k.,  k  ,  etc. 

J  ** 

To  take  proper  account  of  such  correlated  parameters  one  needs  to 
modify  the  Fourier  amplitude  sensitivity  method  analysis.  One  way  of  doing 
this  is  to  incorporate  the  concept  of  correlated  parameters  into  the  formu¬ 
lation  of  the  sensitivity  analysis  theory  from  the  outset.  This  means 
that  one  can  not  use  the  ansatz  (2.3)  according  to  which  the  probability 
density  P(uJ  is  written  as  the  product  of  the  P(u-).  We  have  not  pursued 
this  approach,  and  it  is  not  clear  what  forms  our  theory  will  take  when 
this  simplifying  assumption  does  not  hold.  Research  on  this  problem  is 
clearly  of  interest,  particularly  for  the  application  of  this  method  of 
sensitivity  analysis  to  economic  model  systems. 


Appendix  I 


The  error  e  in  approximating  an  integral  over  e_-space  by  a  line  Integral 
over  the  search  curve,  for  the  Fourier  coefficients  of  frequencies 
Pcc^ ( p=  1 ,2,. . . )  is  represented  by  the  difference 

ipu,„s  ipe.  i  pe  „ 

<f(s)e  ‘  >  -  <f(e)e  £>  =  e(f(o)e  l)  .  (Al.l) 


Here,  e  is  given  by  the  sum 


t  =  l*  cr  5(r-w-p.ojz) 


(A1.2) 


wnere  the  prime  on  the  summation  excludes  the  Cq 
the  r,  ■  r2  -  ...  •  r..,  «  -  ...r„  -  O.r,  ■  p, 

The  Kronecker  delta  is  defined  as 


n  term,  i.e., 

.  •  •  >v 

l 

Fourier  coefficient. 


6  (r-'jj)  = 


1  j  for  r-w  =  0 


0/ 


for  r-u  f  0 


(A1.3) 


This  result  is  obtained  by  using  the  definition,  Eq.  (2.14),  of  f(e)  as  a 
multiple  Fourier  series  in  the  9 ^ 1 s  in  Eq.  (Al.l).  As  can  be  seen  from 
Eq.  (A1.2),  each  time  r-w  =  p  y  one  obtains  a  contribution  to  the  error 
e  made  in  equating  the  s  and  9^  space  integrations.  We  refer  to  the  values 
of  jr  which  satisfy  [r.-t.  =  p  w  as  interferences .  The  weight  of  the  error 
is  just  the  Fourier  coefficient  cr  evaluated  at  the  value  of  r  for  which 
r.*u  *  p  u^.  Since  Fourier  coefficients  tend  to  decrease  in  magnitude  as  their 
index  r  increases,  we  see  that  postponing  the  occurrence  of  interferences 
to  higher  values  of  r  will  lead  to  a  smaller  error  e. 


It  is  convenient  to  define  ,  the  order  of  an  interference  such  that 
r*w  t  p  w  for 

prJ^H.+p,  .  (A1.4) 

i=l  1  1  4 


The  higher  the  order  of  interference,  the  smaller  the  error  e  for  a  given 
output  function. 

The  values  of  £  which  lead  to  interferences  are  controlled  by  the 
choice  of  frequencies  w.  Thus,  a  "judicious"  choice  of  u  leads  to  error 
terms  which  are  small.  By  way  of  illustration  we  return  to  the  two 
dimensional  frequency  choices  given  in  Eqs.  (2.12).  Let  us  say  that  we  are 
interested  in  the  Fourier  coefficient  C ^  and  thus  p£  =  2.  For  the  first 
choice  (<^  =  1 ,  =  2),  the  lowest  interference  arises  when  say  p-j  =  2, 

p 2  =  -1  so,  with  reference  to  tq.  (A1.4),  1.  This  choice  of  frequencies 

leads  to  a  poor  coverage  of  £  space  and  therefore  a  large  error  in  equating 
s  and  £  space  integrals.  The  second  frequency  choice,  =  11  u>2  =  13,  leads 
to  an  interference  when  p-j  =  13  and  p2  =  -11  so  that  =  22.  The  e-space 
coverage  is  much  improved  and  we  obtain  a  more  accurate  value  for  the  e^space 
integration.  By  choosing  frequencies  which  are  increasingly  incommensurate 
as  measured  by  an  increasing  value  of  the  order  of  interference  one  has, 
as  Mj-+  “,  the  strict  equality  of  s  and  e  space  integration. 
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Appendix  II 


The  results  of  Append’ x  I  must  ba  modified  to  account  for  the  use  of 
a  finite  number  of  points  N  that  are  used  to  carry  out  the  s-space  integration. 
The  difference  between  the  s-space  quadrature  and  the  e_-space  average  can  be 
written  as 


where  here  e*  is  found  to  be  (see  ref.  3) 

e*  *  j  1'  cr  6(r-o).-p^j,-jM)  .  (A2.2) 

j = -®  — 

As  in  Eq.  (A1.2),  the  prime  on  the  summation  excludes  the  coefficient 

c  =  c  from  the  summations  over  r.  (i  =  l,2,...,n).  The  error  e*  now 
r  p  i 

—  i 

includes  contributions  when  r.  satisfies 

r-w-P...  =  jiJ  i  j  =  ±1,  ±2 .  (A2.3) 

A,  A 

That  is,  in  addition  to  the  interference  errors  (those  arising  for  j=0), 
each  time  an  integer  multiple  of  N,  the  number  of  points  in  the  quadrature, 
equals  the  frequency  r»u- P^,  an  additional  error,  whose  size  is  given  by 
the  e-space  Fourier  coefficient  of  that  frequency,  also  occurs.  We  term 
these  contributions  to  the  error  as  al iases.  The  aliasing  phenomena,  which 
arise  from  the  finite  number  of  points  used  to  numerically  evaluate  the 
Fourier  coefficients,  can  also  be  controlled  by  appropriate  choice  of  u 
and  N.  Postponing  the  aliasing  to  values  of  cr  with  r_  larger  and  larger  will 
yield  a  more  accurate  Fourier  coefficient  approximant  C*y  .  As  for  inter¬ 
ferences,  we  define  the  order  of  an  alias  such  that  r_-w  f  p  u  +  jN 
(j  *  *1,  ±2,  ..•)  for 
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Postponing  the  occurrence  of  aliasing  by  appropriate  choices  of  u  and  N 
which  lead  to  a  large  value  of  M  decreases  the  error  due  to  aliasing.  For 
further  details  see  ref.  (3). 


(A2.4) 


MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BUREAU  OF  STANDARDS- 1963- A 


*■> 


Appendix  III 


In  the  numerical  computation  of  the  Fourier  coefficients,  an  N  point 
approximation  to  the  t^-space  integration  is  used.  In  order  to  assess  our 
approximations,  we  have  constructed  upper  bounds  on  the  error  made  for  given 
sets  and  classes  of  output  functions  to  be  analyzed.  We  have  described 
this  error  analysis  in  some  detail  in  ref.  (3)  and  will  only  sketch  here  its 
features  and  adopt  them  to  the  new  approach  to  sensitivity  analysis  developed 
in  this  work. 

The  error  made  in  the  numerical  evaluation  of  the  Fourier  coefficients 
is  just  the  difference  between  the  N-point  quadrature  formula  and  the  integral 
over  all  u-space.  For,  if  we  could  use  a  space  filling  search  curve  in  u- 
space,  Weyl 's  theorem  would  be  exact  and  the  Fourier  coefficient  evaluation 
would  be  exact.  Thus,  we  define  the  error  .1 .  as 


A i  ' 


(A3.1 ) 


where  g^te)  =  f(6)e10£,  e  (elq,  e2q . 3nq)  and  ejq  *  Ujsq  with 

j  =  l,2,...n.  We  now  construct  an  upper  bound  A*up  to  the  error  A^  which 
is  a  function  of  u,N  and  the  output  function  to  be  evaluated.  Note 
that  g  (e_)  is  multiply  periodic  in  e_,  and  if  the  partial  derivatives 

X»  " """ 


♦ 

ap  gt(e) 


apl*ap2 
91  02  • 


t 

o<p  £ 


n 

p-1;  o<p  .sp-l;  l 
J  j=l 


(A3. 2) 


are  of  bounded  variation,  then  the  Fourier  coefficients  bp  of  gt(e) 
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defined  by 


9(g.)  *  l  br 

r  - 


r-e 


(A3. 3)  - 


satisfy 

!>,*>,  I  I  <«.4> 

-  V  j  =  1  J 

with  a  independent  of  r  and  rt  =  max(l,jr.|).  The  equality  of  Eq.  (A3. 4) 

r  J  tj 

provides  the  Fourier  coefficients  b^up  for  the  construction  of  asup. 

That  is,  since  g^fej  is  multiply  periodic  in  e_,  we  use  its  multiple  Fourier 
expansion  in  and  obtain 


with 


£  A 


sup 

l 


sup  a 

l  t 


■ftc 

l 

j=-a 


bsup 


6  ( r -oi-jll) , 


(A3. 5) 


(A3. 6) 


One  can  readily  construct  the  function  gsup(ej  whose  Fourier  coefficients 
are  b*up  and,  with  Eq.  (A3.1),  find  £SUp  in  terms  of  oi,N  and  the  properties 
of  the  output  functions  to  b*>  considered.  The  description  of  the  output 
function  is  given  by  the  boundedness  of  its  partial  derivatives  as  described 
by  the  value  of  p*  in  Eq.  (A3. 2).  In  this  fashion  we  have  in  reference  (3) 
constructed  bounds  on  the  error  for  sets  of  w,N  which  we  have  used  in 
numerical  calculations.  The  most  important  conclusion  to  be  drawn  from 
the  error  analysis  presented  in  reference  (3)  is  that  the  Fourier  coefficients 
can  be  calculated  with  good  accuracy  by  judicious  choice  of  N  and  the  frequency 
set  u  and  that  comparison  of  Fourier  coefficients  with  a  "safety  factor" 
of  about  one  order  of  magnitude  yield  valid  sensitivity  measures. 
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The  bounds  on  the  Fourier  coefficients  can  readily  be  combined  to  yield 

bounds  on  the  partial  variances  S  .  Thus,  the  Fourier  coefficient  analysis 

“a 

Is  directly  applicable  to  the  partial  variances.  It  is  Important  to  note 
that  as  one  examines  the  Fourier  coefficients  of  higher  harmonics  of  a 
fundamental,  the  accuracy  with  which  they  are  approximated  degenerates 
since  one  is  using  the  same  number  and  placement  of  points  In  a  quadrature 
formula  Involving  a  function  which  is  increasingly  more  oscillatory. 

However,  the  s-space  Fourier  coefficients  themselves  fall  off  In  magnitude 
as  one  works  with  higher  frequencies  so  that  the  Fourier  coefficients  for 
the  higher  frequencies  do  not  have  to  be  as  precisely  approximated  as  those 
for  the  lower  frequencies. 


Appendix  IV 


2 

The  part  of  the  total  variance  o  arising  from  the  uncertainty  of  the 
f'th  parameter,  when  the  out^t  function  is  averaged  over  the  uncertainties 
in  all  other  parameters,  is  found  by  first  integrating  over  all  parameter 
uncertainties  except  for  the  n'th.  This  integration  is  best  done  in  space 
where  one  may  use  the  myltiple  Fourier  decomposition  of  Eq.  (2.14).  Thus 


f(<g  =  /de1...de1-1de1+1...denf(e1,...en) 


l  C00...pr..00ex^1Pi6^- 


(A4.1) 


Now  form  the  variance  o,  of  f(e  )  with  respect  to  the  uncertainty  in  k 


°\  *  /f2(0i)d9a  ’  J2 


(A4.2) 


n  l  JcOO...p0...Ool 

p  *-*  rl 


We  have  found  before  (cf.  Eq.  (2.18))  that  the  e-space  Fourier  coefficient 

o 

c00  p  00  eC|ua*s  s -space  coefficient  so  that  can  also  be 
written  as 


2 

o. 


(A4.3) 


As  before,  we  must  modify  this  result  to  correspond  to  the  finite  summations. 
We  define  <?*  to  be 


p=-(N/2-l ) 


(A*  + 

P«, 
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2* 

Now  form  the  quotient  arising  from  the  i'th  parameter  and  the  total 
*2  * 

variance  o  and  call  it  S  ,  the  partial  variance 

“i - 


N/2 


*  2 

I C  1 

5*  -  P=-(N/2-l)  p^' 

V  '  N/2  *  2 


(A4. 


j=-(N 


1/2-1)' 


1C, 


This  result  is  Eq.  (2.25)  of  the  text. 


Appendix  V 


We  wish  to  establish  the  symmetry  relations 

fU/2  +  s)  =  fU/2  -  s)  (A5.1 ) 

f(-*/2  -  s)  =  f (-ir/2  +  s)  (A5.2) 

of  the  output  function  f(s)  for  -ir/2  i  s  s  ir/2  and  to  see  how  these 
relations  reduce  the  size  of  the  discrete  samples  along  the  closed  search 
curve. 

To  do  so,  we  stipulate  that  the  .^'s  are  odd  integers  and  set  *  2k  +  1 

Thus 

sin[(2k+l )]  =  sin[(2k+l  ){tt-s)]  =  -  sin[(2k+l  )(n+s)]  =  -  sin[(2k+l ) (2-ir-s ) ] . 

(A5.3) 

Applying  these  identities  to  all  the  parameters,  we  establish 


f(ir-s)  =  f(sinto-|  (tr-s),  sin^  (ir-s),...) 

*  f(sinw-|S  ,sinw2s ,. . . )  *  f(s) 

f (tt+s )  =  f(-sin  ^s.-sinu^s,...)  =  f(2ir-s)  . 


(A5.4) 

(A5.5) 


Since  f(s)  *  f(s+2ir)  by  construction,  eq.  (A5.5)  also  implies 

f  ( -TT+s )  =  f  ( s )  (A5.6) 

Equation  (A5.4)  states  that  f(s)  in  the  quadrant  tt/2sks:it  is  the  mirror 
image  about  ir/2  of  f(s)  in  the  quadrant  Osssir/2;  eq.  (A5.6)  states  that  f(s) 
in  the  quadrant  -it£s*t/2  Is. the  mirror  image  about  -ir/2  of  f(s)  in  the 
quadrant  -ir/2sss0. 

Thus  we  need  only  evaluate  f(s)  in  the  range  =ir/2if (s)sn/2.  This  reduces 
the  computational  load  by  a  factor  of  two  relative  to  what  was  stated  in 
references  (1)  through  (3). 
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Appendix  VI 


The  working  equations  of  our  sensitivity  method  yield  the  Fourier 
coefficients  of  Eqs.  (3.13).  Here  we  show  how  to  obtain  these  equations  by 
use  of  the  symmetries  of  the  output  functions  f(s)  in  s,  Eqs.  (3.8),  and 
the  symmetries  of  the  trigonometric  functions. 

The  Fourier  coefficients 


Aj  ‘  r  X  fq  C0S  jSq 


B3  '  F  qI,  fq  S1"  3Sq 


(A6.1) 


can  be  written,  using  these  symmetries,  as  sums  over  the  half  Interval 
n/2  i.  sq  <  tt/2 

AJ '  ^  q-(Li)/2tcosjsq  *  “sjs-r-q]f< 

,  (r-1 )/2 

+  I  l  [cosjs  +  cosjs  ]f  (A6.2) 

q=l  M  H  H 

1  0 

i  I  Csinjs  +  sinjs_  ]f 

J  r  q*-(r-1 )/2  q  rq  q 

,  (r+1 )/2 

+  '7  q*i  Cs1njsq  +  SinjSr-q]fq  (A6-3) 

where  we  now  define  sq  *  wq /r  and  set  f(Sq)  *  fq. 

Application  of  the  addition  theorems  for  trigonometric  functions  allows 
us  to  rewrite  eqs.  (AC. 2)  and  (A6.3)  as 


(A6.3) 


i  <  (r-l)/2 

Aj.l[i  M-D%0  ♦  l}  Cfjf.jlcMaa 


(A6.4) 
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i  <  (r-1 )/2 

Bf  =  kl-(-Dj]  l  [frf  J  sin^  . 

J  q=l  J  "« 


(A6.5) 


Recalling  that  r  *  2q+l  we  obtain  the  equations  (3.13)  given  in  section  3. 
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Figure  1:  Definition  of  Terms  for  Sirulation  Models. 


Figure  2 


Figure  3 


Figure  4 


Figure  5 


Figure  6 


Figure  7 


Figure  8 


9,-space  coverage  with  frequencies  W|  =  1 ,  *  2  and  100  , 

quadrature  points  N.  The  lines  are  obtained  from  e  *  «  s(mod  2n) 

X*  X 

0  s  s  <  2it  and  the  points  from  0„  =  u  s  (mod  2w)  with 

tq  i  q 

Sq  a  2wq/N;  q  =  1,2,...N. 

9-space  coverage  with  =  11 ,  =  I3  an< i  N  =  100  computed 

in  the  same  fashion  as  in  Figure  1. 

Zero  power  gains  versus  tine  for  the  v  =  2-+1 ,  v  *  3+2,  and  v  •  4+3 
vibrational  transitions  of  HF.  Curves  show  the  time  histories  for 
both  the  nominal  parameter  values  and  for  the  statistical  mean 
value  (average  over  the  parameter  uncertainties). 

Coefficients  of  variation  of  the  gains  on  v  *  2+1 ,  v  s  3+2,  and 
v  *  4+3  vibrational  transitions  of  HF,  as  functions  of  time. 
Partial  variances  for  the  zero  power  gain  of  the  v  =  2-+1  band 
of  HF.  Curves  are  labeled  by  the  rate  coefficients  to  which 
they  correspond;  [F]g  curve  is  partial  variance  due  to  uncer¬ 
tainty  in  initial  F  atom  concentration. 

Partial  variances  for  the  zero  power  gain  of  the  v  *  3+2  band 
of  HF.  Curves  are  labeled  by  the  rate  coefficients  to  which 
they  correspond;  [F]g  curve  is  partial  variance  due  to  uncer¬ 
tainty  in  initial  F  atom  concentration. 

Partial  variances  for  the  zero  power  gain  of  the  v  *  4+3  band  of 
HF.  Curves  are  labeled  by  the  rate  coefficient  to  which  they 
correspond;  [F]g  curve  is  partial  variance  due  to  uncertainty  In 


initial  F  atom  concentration. 


Figure  9:  World  population  versus  time. 

Figure  10:  Partial  variances  for  world  population.  Curves  are  labeled  by 


parameters  to  which  they  correspond.  Three  parameters,  POLN,  BETA 
and  C  have  neglibible  partial  variances  prior  to  2020.  After 
that  time  they  have  partial  variances  about  equal  to  those  for 
BRN  and  DRN.  They  are  omitted  from  the  figure  so  as  not  to 
complicate  it.  Other  parameters  have  negligible  partial  variances. 

Figure  11:  Quality  of  life  versus  time. 

Figure  12:  Partial  variances  of  quality  of  life.  Curves  are  labeled  by  the 
parameters  to  which  they  correspond.  The  parameters  POLM,  BETA, 
arH  NR I  have  partial  variances  similar  to  CIDN.  These  are  not 
shown  in  the  figure  to  avoid  confusion.  Other  parameters  have 
negligible  partial  variances. 


-104- 


. 


CO 

Ul 

_l 

03 

to 

< 

LU 

►4 

3 

CC 

_J 

< 

< 

> 

> 

1- 

LU 

z 

O 

LU 

a 

Ul 

z 

L5 

Ul 

Z 

a. 

< 

LU 

CC 

a 

z 

to 

in  a 
lu  to 

z  z 

LLt  O 

>  a  ~ 
a  —  to  H- 
H  H-  H  U 
3  O  J  « 
0-  LU  3  C3 

I—  U_  L0  LU 

3  LL  uj  o; 
O  UJ  CC  Cc 


H 

CO 

z 

UJ 

LU 

_l 

a 

03 

z 

< 

LU 

1— 1 

CL 

QC 

LU 

< 

a 

> 

to 

ce  to 

LU  U1 
H-  3 
UJ  -I 

s:  < 
<  > 
0 i 

<  u. 
CL  O 

-J  Ul 

o  o 
cc  z 
t-  < 
z  cr 

O  W 

u 


to 

> 

0£ 

t- 

Ul 

z 

►“ 

44 

LU 

< 

s: 

< 

QC 

cc 

Ul 

< 

o 

Q. 

z 

3 

< 

-J 

1- 

< 

z 

u 

LU 

£ 

H* 

44 

to 

Q£ 

44 

Ul 

t- 

CL 

< 

X 

h 

LU 

to 

mma 

CO 

> 

J— 

1- 

z 

z 

Ul 

44 

C 

u 

h- 

cc 

U. 

LU 

LL 

t_> 

LU 

z 

o 

3 

o 

-1 

z 

< 

o 

a 

44 

44 

to 

h 

to 

to 

Ul 

44 

cc 

1- 

a 

< 

Ul 

J— 

o: 

to 

_ 

W 

03  Z.  > 

LU  t  U- 

£  5  ° 

z  §  UJ 

—  go 

U.  u  Z 

UJ  < 

a  cc 

Z  W 


to  LU 
UJ  X 
-J  t- 
03 

<  o 
—  z 
cc  — 

<  cc 

>  =3 

a 

t— 

Z  LU 
LU  O 

a  z 
z  < 
uj  x 
a.  c_> 

LU 

a 


to  _J 
Ul 

<  c 
o 

u-  x 

o 

LU 
Ul  X 
to  J— 
cc 

3  U. 

o  o 

O 

z 

LU  3 
X  cc 
H 

o 

tr  t- 

z 

—  z 
a:  3 
3  CC 

a 

X 

ui  o 
o  ce 
z  u. 

< 

X  Ul 

u  o 
••  z 
to  I-  < 
cc  o  x 

LU  Z  O 
H 

Ul  O  O 
E  C  C 
< 
cc 
< 
c_ 


^ 


Ho  NOT  CHANGE  DURING  THE  COURSE  OF  A  SINGLE  RUN  OF  THE  MODEL. 
Do  NOT  CHANGE  FROM  RUN  TO  RUN  OF  THE  MODEL. 
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